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ABSTRACT 

This  thesis  concerns  developing  a  program  that  takes  an 
aerial  photograph,  and  a  set  of  Digital  Terrain  Elevation 
Data  (DTED)  that  is  defined  over  the  area  of  the  photograph, 
and  generates  a  synthesized  view  that  represents  what  a 
camera  would  see  from  a  different  location.   The  elevation 
data  points  are  grouped  into  triangular  panels  that  are 
projected  to  the  reference  image  by  three  dimensional 
transformation  equations.   Shading  for  the  synthesized  image 
is  determined  from  the  reference  image.   The  pixels  of  the 
reference  image  that  fall  within  a  triangular  panel  are 
collected  and  averaged.   When  a  new  observer  location  is 
selected,  the  panels  are  projected  to  the  new  synthesized 
image  plane.   A  z-buffer  approach  and  a  polygon  fill 
algorithm  were  used  to  remove  hidden  surfaces  of  the 
synthesized  view. 

This  program  is  tested  on  both  artificial  and  real  data. 
Dther  characteristics  and  performance  measurements  of  the 
program  are  also  analyzed  here.   The  quality  of  the  syn- 
thesized image  from  real  data  was  affected  by  the  low 
resolution  of  the  terrain  elevation  data,  and  yielded  less 
desirable  results  than  could  be  expected  of  a  higher 
resolution  terrain  model  . 
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I .   INTRODUCTION 

A.   COMPUTER  IMAGE  GENERATION  FROM  AERIAL  PHOTOGRAPHY 
The  main  objective  of  this  study  was  to  develop  a 
program  that  takes  a  digital  photographic  image  and  a  file 
of  terrain  elevation  points  defined  over  that  image  as 
input,  then  produces  as  an  output  a  synthesized  perspective 
view  .   The  synthesized  view  is  a  rotated  3-dimensional  C3D) 
perspective  representation  of  the  original  photographic 
image.   The  main  application  of  this  study  is  to  generate  a 
different  perspective  of  a  terrain  model .   This  may  be  used 
to  generate  different  views  that  a  pilot  of  an  aircraft 
could  expect  following  different  flight  paths  through  the 
same  area.   Further  study  may  make  it  feasible  to  generate 
synthesized  images  fast  enough  to  simulate  a  real  time  image 
display  of  a  flight  for  a  mission  briefing  or  to  be  used  as 
a  training  aid.   Another  application  could  be  for  training 
men  on  optically  guided  missiles.   With  high  resolution 
images,  a  flight  path  through  a  battlefield  could  be 
simulated  that  would  have  all  the  visual  characteristics  of 
an  actual  flight  without  the  expenditure  of  a  live  missile. 
The  generation  of  a  shaded  image  as  a  3D  picture  provides 
unique  problems  for  3D  graphic  displays.   The  data  which 
comprises  a  photographic  image  consists  of  an  array  of 
pixels,  each  of  which  has  a  defined  grey  level  or  shade. 
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There  are  ESS    different  levels  of  grey  that  may  be  assigned 
to  a  pixel .   This  study  concerns  taking  a  photographic 
perspective  image  and  a  S-dimensional  (ED)  array  of  eleva- 
tions defined  in  a  grid  covering  the  area  of  the  image  as 
inputs.   The  grid  of  elevations,  called  the  terrain  model, 
is  geometrically  related  to  the  photographic  image  through  a 
perspective  projection  transformation  that  equates  the  world 
coordinates  of  the  elevation  points  to  the  object  coor- 
dinates of  the  image.   A  synthesized  image  from  a  different 
observer  location  is  then  generated.   The  new  synthesized 
view  should  approximate  what  would  be  seen  by  a  camera  from 
the  new  observer  position. 

The  differences  between  the  original  and  synthesized 
images  will  be  affected  by  the  resolution  of  both   the 
photographic  image  and  the  terrain  models.   Higher  resolu- 
tion of  the  original  models  will  result  in  a  closer  approx- 
imation in  the  synthesized  image.   Another  complication  or 
ambiguity  arises  when  details  which  should  show  up  in  the 
synthesized  view  were  not  present  in  the  original  image.   A 
method  must  be  devised  so  that  it  can  fill  in  areas  which 
became  visible  in  the  synthesized  image  that  were  hidden  in 
the  original  reference  image.   The  solution  to  this  hidden 
surface  problem  is  further  addressed  in  the  discussion  of 
the  grey  scale  referencing  algorithm . CRef .  13 
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B.   TERRAIN  ELEUATION  AND  PHOTOGRAPH  AS  INPUTS 

In  this  study  a  high  altitude  aerial  image  of  MoFfett 
Field,  California  was  used  as  the  original  reference  photo- 
graph.  The  photographic  image  uas  supplied  by  the  Defense 
flapping  Agency  (DMA)  and  had  a  resolution  of  approximately  1 
meter  per  pixel .   The  terrain  model  corresponding  to  the 
reference  image  uas  provided  as  Digital  Terrain  Elevation 
Data  (DTED)  by  the  DMA  and  consisted  of  elevation  points 
taken  every  second  of  a  degree  change  in  latitude  and 
longitude.   This  gives  an  approximate  resolution  of  30 
meters  per  elevation  point  in  a  north-south  direction  and  23 
meters  per  elevation  point  in  an  east-west  direction. 

The  synthesized  view  was  restricted  to  a  northerly 
direction  which  simulates  an  aircraft  flying  from  south  to 
north  with  the  image  plane  perpendicular  to  the  direction  of 
flight.   To  allow  for  different  flight  patterns  would 
require  development  of  an  algorithm  that  would  provide  for 
rotation  of  the  image  coordinates  which  is  beyond  the  scope 
of  this  study  .   The  main  idea  is  to  generate  a  synthesized 
view  that  is  rotated  from  the  original  photograph  by 
approximately  SO*  and  explore  the  concepts  of  the  algorithms 
required  to  da  this.   Although  speed  was  not  a  major  issue, 
the  size  of  the  terrain  model  was  limited  to  a  50  x  50  grid 
array,  or  2500  data  points,  to  minimize  the  time  for 
synthesized  image  generation. 
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C.   ALGORITHM  ISSUES 

1 •   Greu  Scale  Referencing 

To  determine  the  grey  scale  values  of  the  pixels 
that  make  up  cur  synthesized  visu,  the  terrain  model  is 
first  divided  into  triangular  panels.   The  vertices  of  the 
triangular  panels  are  then  mapped  into  the  original  DMA 
image  using  the  perspective  projection  transformation  that 
projects  georectangular  coordinates  into  reference  image 
coordinates.   The  pixel  grey  scale  values  that  fall  within 
the  projected  triangular  panel  are  then  averaged.   The 
average  grey  scale  value  is  permanently  assigned  to  that 
particular  panel  .   Uhen  the  synthesized  image  is  constructed 
the  triangular  panel  is  mapped  to  the  new  image  vieui  and 
filled  uith  the  assigned  average  grey  scale  value.   In  this 
way  the  sample  of  pixels  that  fall  within  the  triangular 
panel  are  mapped  from  the  original  reference  image  to  the 
synthesized  image.   This  method  of  mapping  the  triangular 
panels  to  the  synthesized  image  also  salves  the  problem  of 
assigning  a  grey  scale  value  to  hidden  surfaces  of  the 
reference  image  because  they  are  automatically  assigned  the 
value  of  the  surrounding  pixels.   Since  the  average  grey 
scale  value  for  a  triangular  panel  is  dependant  upon  the 
resolution  of  the  terrain  model ,  the  grey  scale  value 
assigned  to  the  hidden  surface  will  also  be  affected 
CRef.  13. 
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The  smaller  the  triangular  panels  are,  the  smaller 
the  area  that  must  be  collected  and  averaged  in  the  original 
image.   This  means  a  much  better  synthesized  view  can  be 
constructed  that  contains  more  of  the  attributes  of  the 
original  image.   For  this  reason  the  resolution  of  the 
terrain  and  reference  image  model  is  very  critical  to 
obtaining  an  accurate  synthesized  view.   Using  the  resolu- 
tion of  the  terrain  and  image  models  used  in  this  study,  the 
approximate  number  of  pixels  that  must  be  averaged  in  the 
reference  image  for  each  triangular  plane  would  be  1/5  (30m 
x  23m .1(1  pixel/m)  ■  345  pixels.   This  is  very  coarse  and 
does  not  allow  for  optimal  generation  of  the  synthesized 
view  . 

2  •   Hidden  Surface  Elimination 

There  are  surfaces  that  may  be  discernible  in  the 
reference  image  but  become  hidden  in  the  synthesized  view. 
The  z-buffer  algorithm  was  used  to  accomplish  the  hidden 
surface  elimination.    The  z-buffer  is  an  array  that 
contains  the  depth  or  distance  to  the  observer  location  for 
each  pixel  that  is  to  be  visible  in  the  synthesized  image. 
As  each  triangular  plane  is  mapped  to  the  synthesized  view 
the  location  and  depth  of  each  pixel  within  the  plane  is 
determined.   The  depth  of  the  pixel  to  be  written  at  a 
certain  location  is  compared  to  the  depth  of  any  pixel  that 
may  have  been  previously  written  to  the  same  location.   If 
the  depth  or  distance  of  the  new  pixel  to  the  observation 


14 


point  is  shorter  than  the  previous  pixel,  its  depth  is 
written  to  the  z-buffer  and  the  grey  scale  value  oF  the 
pixel  is  placed  into  the  synthesized  image.   IF  the  depth  is 
larger,  no  updating  occurs  and  a  new  pixel  is  obtained  in 
the  process.   The  z-buFFer  works  very  closely  with  the  next 
algorithm  to  be  considered.   CReF.  2,  pp.  265-2672 
3 .   Polugon  Fill  Algorithm 

Screen  coordinates  are  generated  For  the  three 
vertices  oF  each  triangular  panel  as  it  is  mapped  to  the 
synthesized  image.   Screen  coordinates  are  designated  as  IA 
and  JA  values  with  the  IA  values  representing  the  columns 
and  the  JA  values  the  rows.   The  location,  IA,JACO,0}, 
designates  the  upper  leFt  hand  corner  oF  the  screen  and  the 
maximum  screen  coordinate,  I A , JAC51E, 512} ,  the  lower  right 
hand  corner.   An  active  edge  list  CAEL}  is  generated  by 
computing  a  line  between  each  oF  the  translated  vertices. 
For  each  line,  the  IA  coordinate  corresponding  to  the 
maximum  JA  value  oF  the  line,  the  amount  IA  changes  For  each 
one  unit  step  oF  JA,  and  the  total  span  oF  JA  are  stored 
into  the  AEL .   The  three  lines  generated  For  each  translated 
triangular  plane  will  Form  another  closed  triangle.   By 
using  the  parameters  stared  in  the  AEL  the  location  oF 
pixels  enclosed  by  the  translated  triangular  plane  can  be 
determined,  and  the  corresponding  array  points  within  a 
Frame  buFFer  are  changed  From  a  0  to  a  1 . 
CReF.  2,  pp.  76-733 
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After  all  of  the  enclosed  pixels  have  been  marked 
within  the  frame  buffer,  the  buffer  is  scanned  row  by  raw. 
If  a  value  of  1  is  found,  then  the  depth  is  calculated  for 
that  point  and  compared  with  the  depth  value  stored  in  the 
z-buf f er .   If  the  depth  value  is  smaller,  that  pixel  is 
located  closer  to  the  observer  location  and  the  grey  scale 
value  for  that  pixel  is  written  to  the  synthesized  image 
file.   As  can  be  seen  the  fill  and  hidden  surface  algorithms 
work  together  to  generate  the  new  image.   The  implementation 
of  these  algorithms  are  explained  in  further  detail  later. 

In  Chapter  II  there  will  be  a  discussion  of  basic 
photographic  geometry  to  develop  an  understanding  of  the 
transformation  equations  necessary  to  map  object  coordinates 
into  image  coordinates  and  for  image  plane  rotation. 
Chapter  III  will  detail  program  considerations  based  on 
image  and  elevation  data  as  well  as  the  algorithms  used  to 
generate  the  synthesized  view.   Chapter  IU  will  discuss 
possible  ways  to  improve  the  transformation  program   and 
discussion  of  topics  for  passible  further  study.   An  outline 
of  the  program  is  contained  in  Appendix  A  that  gives  a  short 
discussion  of  each  subroutine  as  to  its  purpose,  input  and 
output,  modules  that  called,  and  modules  that  reference  the 
subroutine.   Appendix  B  contains  the  entire  3D  transforma- 
tion program. 
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I  I  •   PHOTOGRAPHIC  GEOMETRY 

A.   BACKGROUND 

To  understand  many  of  the  concepts  used  in  this  study,  a 
basic  background  in  photographic  geometry  is  presented.   The 
relationship  between  the  image  space  and  abject  space  is  the 
basis  For  many  of  the  equations  that  help  to  generate  the 
reference  and  synthesized  images  for  visual  display  .   The 
objective  of  this  chapter  is  to  present  the  concepts  that 
allow  the  transformation  of  3D  objects  to  a  ED  image  and  the 
parameters  evolved. 

1  .   Perspective  and  Parallel  Projection 

A  parallel  projection  is  a  projection  in  which  the 
projection  lines  from  the  abject  to  the  image  plane  never 
converge.   When  an  object  is  viewed  by  parallel  projection, 
its  size  would  never  change  as  the  camera  is  moved  closer  or 
further  away.   In  contrast,  a  perspective  projection  has  all 
the  projection  lines  from  an  object  converge  to  a  perspec- 
tive center.   A  perspective  projection  imitates  how  we  see 
things.   An  example  would  be  a  picture  of  railroad  tracks. 
The  tracks  would  appear  to  became  closer  together  when 
further  away  from  the  observation  point.   In  a  parallel 
projection  the  tracks  wculd  be  the  same  distance  apart  along 
the  entire  length.  CRef.  3,  pp.  133-134] 

Since  a  camera  is  generally  designed  to  photograph  a 
rather  large  area,  it  involves  perspective  projection.   The 
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camera  view  represents  what  an  observer  would  see  standing 
at  the  same  location,  and  the  images  generated  are  perspec- 
tive images.   This  means  that  the  equations  used  to  trans- 
form the  object  space  into  the  image  space  must  be  perspec- 
tive transformations. 

2 .   Image  Coordinate  Space 

The  image  plane  is  the  plane  of  the  photograph  to 
which  the  object  points  are  mapped.   It  has  a  ED  coordinate 
system  to  which  each  point  of  a  3D  object  is  translated  to 
(Fig.  2.1).   The  indicated  principal  point  CIPP)  is  the 
center  of  the  image  plane  and  has  the  coordinates  of 
Cx,y,0).   The  x,  y,  and  z  axis  represent  a  right  handed 
plane  and  the  perspective  center  CD  lies  along  a  line 
parallel  to   the  z-axis;  then  a  perpendicular  line  is  drawn 
from  L  to  the  image  plane.   The  point  at  which  this  perpen- 
dicular line  intersects  the  image  plane  is  called  the 
principal  point  Co).   This  offset  of  the  principal  point 
from  the  IPP  is  compensated  for  in  the  transformation 
equations  by  xo  and  yo .   The  focal  length  of  the  plane  is 
defined  as  the  distance  from  the  principal  point  to  the 
perspective  center.   For  the  image  plane  coordinate  system, 
each  object  point  CAD  is  graphed  to  a  corresponding  image 
plane  point  CaD  located  at  Cxa,ya,OD,  and  the  perspective 
center  or  focal  point  is  located  at  Cxo,yo,fD.  CRef.  3,  pp. 
135-136] 
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In  generating  a  synthesized  view,  the  perspective 
center  may  be  placed  at  any  location  desired  with  reference 
to  the  image  plane.   It  is  therefore  desirable  to  select  a 
point  along  the  z-axis  such  that  xo  and  yo  become  0.   This 
will  decrease  the  number  of  calculations  required  in 
generating  the  synthesized  view. 
3 .   Dbiect  Coordinate  Space 

The  observer  location  and  each  object  point  position 
is  described  in  world  coordinates,  called  georectangular 
coordinates,  of  X,  Y,  and  2.   The  center  of  the  earth  is 
given  as  (0,0,0),  the  Z-axis  points  directly  to  true  north, 
the  X-axis  points  to  the  intersection  of  0*  Latitude  and  0* 
Longitude,  and  the  Y-axis  the  intersection  of  0*  Latitude 
and  90"  E.  Longitude.   The  principal  or  focal  point  that 
would  describe  the  observer  location  in  georectangular 
coordinates  is  XL,  YL,  and  ZL .   Each  object  point  CA) 
located  in  the  object  space  is  identified  by  XA ,  YA,  and  ZA 
as  shown  in  Figure  2.2.   CRef  .  3,  p.  135D 

B.   IMAGE  PLANE  ROTATION 

To  align  the  x,  y,  z  coordinates  of  the  image  plane  to 
the  desired  viewing  direction  requires  rotation  about  the  X, 
Y,  and  Z  axis  of  the  georectangular  coordinate  system.   In 
general  to  transform  one  3D  coordinate  system  requires  a 
matrix  multiplication  of  the  form  A  -  CHUB.   The  A  repre- 
sents a  vector  in  the  image  space  with  x,  y,  and  z  coor- 
dinates, and  B  is  a  vector  in  the  georectangular  coordinate 
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Fig.  £.1   Image  Plane  Coordinates  CRef.  3,  p.  135) 


Fig.  2.2   Object  Plane  Coordinates  CRef.  3,  p.  136) 
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system  with   X,  Y,  and  Z  components.   This  may  be  written  as 


where 


n  - 


a   b   c 
d   e   f 

g   h   1 


a   b   c 
d   e   F 

g   h   1 


X 

Y 
Z 


ca.n 


CosxX  CosxY  CosxZ 
CosyX  CosyY  CosyZ 
CoszX  CoszY  CoszZ 


C2.2? 


This  maps  the  vector  X,  Y,  Z  to  the  image  space  x,  y,  z. 

The  M  matrix  is  derived  From  the  deFinition  oF  the  direction 

oF  a  vector  A  given  by 


A 


AAA 

Coso  1  +  CosB  j  +  Cost  k 


Ax 
IA 


A 

1   + 


Ay 


fiz   A 

A  A 

1  k 


CE.3D 


A     A  A 

where  i,  j,  and  k  are  unit  vectors  oF  the  particular 
coordinate  system  in  which  vector  A  is  contained.   The 
quantities  ex,  B,  t  are  the  angles  that  vector  A  makes  with 
the  x,  y,  and  z-axis  respectively.   Since  there  are  three 
vector  components  in  the  X,  Y,  Z  system,  each  one  must  make 
its  own  transformation  into  x,  y,  and  z  and  thereby  Farming 
the  3x3  matrix  oF  M.   To  translate  x,  y,  z  into  X,  Y,  Z 
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the  inverse  of  the  M  matrix  is  taken  giving  B  =  CHD'A. 
CRef.  3,  p.  1333 

If  we  can  define  three  orthogonal  vectors  in  georect- 
angular  coordinates  as  R,  S,  and  T  that  mould  describe  the 
desired  viewing  position  of  our  image  plane,  the  M  matrix  is 
easily  derived  by 


Sx 
S| 

Rx 

Rl 
Tx 

T| 


Sy 

S! 
Ry 

Rl 

Ty 
Tl 


Sz 

5I 
Rz 

Rl 
Tz 

Tl 


C2.4) 


Generally  the  image  plane  rotation  is  expressed  in  omega 
CuO,  phi  (3D,  and  kappa  (k).   The  w  is  the  rotation  about 
the  X-axis,  the  JO  is  rotation  about  the  Y-axis,  and  k  is 
about  the  Z-axis.   If  we  rotated  first  about  the  X-axis  by 
an  angle  of  w  radians,  the  terms  of  the  M  matrix  would 
equate  as  fallows 

CasxX  =  CosCO')  -  1 

CosyY  =  CasCw) 

CosyZ  -  CasC90#-wD  -  SinCwD 

CoszY  =  CosCw  +  30')  -  -SinCw) 

CaszZ  ■  Cos(w) 
All  other  terms  equate  to  cos  30*  -  0,  therefore  the  M 
matrix  becomes 
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n  = 


10      o 

0    CosCu)   SinCwJ 
0   -SinCw)   CasCw) 


C  2  .  5  D 


Similarly  a  rotation  32  about  the  Y-axis  produces 


n  = 


CosCffiD 

0 

0 

1 

SinCfflD 

0 

-SinCffl) 

0 

Cos  CD) 


C2.6) 


and  For  a  rotation  k  about  the  Z-axis 


n  = 


CosCk) 
-SinCk) 
0 


SinCk)     0 
CosCk)     0 

0        1 


C2.7D 


By  multiplying  all  three  matrices  together  we  derive  the 
overall  M  transform  in  w,  0,  and  k, 


n    = 


CosCflnCosCkO       CosCw^SinCk)    +    SinCui)Sin(ffl)CosCk) 
-CosCanSmCk)       CosCwDCosCk)    -    SinCwDSinC  ID  DSinCk) 
SinCffl)  -    SinCwDCosCan 


SinCw^SinCk) 
SinC w JCosCkD 


CosCu)Sin(I)CosCk) 
CosCuj)SinCB)Sin(k) 
CosC w)Cos(  0  3 

C2.BD 


This  is  the  general  form  of  the  matrix  that  maps  the 
georectangular  coordinates  to  the  image  plane.   CRef  .  3,  pp 
597-600J 
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The  general  farm  af  the  transformation  matrix  is  used  to 
initially  map  the  elevation  terrain  model  into  the  original 
reference  image.   The  w,  CD,  and  k  mere  supplied  with  the 
original  DMA  photograph  and  represents  the  rotation  of  the 
reference  image  plane  with  respect  to  the  georectangular 
coordinates  at  the  time  the  picture  was  taken.   When  we 
generate  the  synthesized  view  the  image  plane  must  be 
rotated  to  the  desired  viewing  angle  of  the  observer 
(northerly  direction  in  particular).   This  means  that  a  new 
w,  3),  and  k  must  be  calculated,  or  by  defining  new  image 
plane  coordinate  axes  in  terms  of  the  georectangular 
coordinates,  we  can  calculate  the  terms  of  the  M  matrix 
directly  using  the  preceding  equations.   This  is  discussed 
further  in  the  next  chapter. 
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III.   ALGORITHM  CONSIDERATIONS 

In  this  chapter  an  in-depth  analysis  of  the  original 
image  and  terrain  data  is  presented.   This  includes  how  the 
data  was  referenced,  the  size  of  the  data  Files,  how  the 
data  was  used,  and  how  the  elevation  and  image  data  compared 
with  one  another.   The  process  of  translating  from  the 
object  space  coordinates  to  image  space  coordinates  and  then 
to  screen  coordinates  is  considered,  and  the  equations  used 
are  given . 

The  referencing,  fill,  and  z-buffer  algorithms  also  are 
discussed  in  detail.   How  the  data  generated  from  these 
algorithms  is  used  and  put  together  to  produce  the  syn- 
thesized view  will  be  presented.   Any  problems  that  were 
encountered  and  the  eventual  solutions  will  be  discussed  in 
the  appropriate  section  to  which  they  pertain. 

A.   REFERENCE  IMAGE  DATA 

The  picture  of  Moffett  Field  supplied  by  the  Defense 
Mapping  Agency  CDMA},  was  49S3  by  4397  pixels  in  size  and 
came  with  both  a  left  and  right  image.  Only  the  left  image 
was  used  to  generate  the  perspective  views  in  this  study . 
Each  individual  pixel  within  the  original  image  is  desig- 
nated by  coordinates  I  and  J,  where  I  is  the  pixel  column 
and  J  is  the  scan  row.   The  geographic  northeast  corner  has 
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the  I,  J  coordinates  of  CO, 03,  and  the  southwest  corner 
C4997,4999) . 

The  image  display   devices  used  were  capable  of  dis- 
playing images  that  were  only  51E  by  512  pixels  in  size, 
therefore,  the  original  image  was  divided  into  appropriate 
blacks  suitable  for  viewing  called  Frames.  Each  Frame 
contains  an  image  that  is  512  by  512  pixels.   The  coor- 
dinates oF  each  Frame  has  a  Four  digit  I_Frame  and 
J_Frame  value,  and  is  Further  identiFied  as  a  leFt  or  right 
CL  or  R)  image.   The  I  Frame  and  J_Frame  coordinates  are 
designated  in  multiples  oF  512   which  deFine  the  column  and 
row  location  oF  each  Frame.   There  are  10  Frame  columns  and 
10  Frame  rows  with  assigned  coordinates  oF  0000  through 
450B .   To  identiFy  a  particular  Frame  oF  the  Drift  image  one 
First  designates  whether  it  is  a  LeFt  or  Right  image  and 
then  give  the  I  _Frame  and  J  Frame  coordinates.   As  an 
example  L05121024  would  designate  a  LeFt  image  from  the 
second  column  and  third  row.   The  First  Frame  L00000000 
starts  in  the  southeast  corner  and  the  last  Frame  L4B0B450B 
is  in  the  northwest  corner.   The  disparity  between  the 
starting  location  oF  the  Frame  coordinates  and  the  I  and  J 
coordinates  oF  the  original  image  must  be  compensated  For  in 
the  equations  that  are  used  to  determine  the  location  oF 
individual  pixels  within  a  Frame  image  as  shown  later. 
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1  .   Object  to  Reference  Image  Transf ormatian 

Every  object  point  is  converted  From  its  3D  X,  Y, 
and  2  georectangular  coordinates  into  the  ED  x  and  y  image 
coordinates  using  the  A  ■  CnUB  equation  discussed  earlier. 
The  vector  from  the  perspective  center  to  each  object  point 
is  defined  by  CXA-XL),  CYA-YL),  and  CZA-ZLD.   This  vector  is 
mapped  into  the  reference  image  plane  coordinates  of  x ,  y, 
and  2.   Since  every  vector  in  the  image  plane  is  directed 
from  the  perspective  center  or  focal  point  to  the  image 
point  (a),  the  z  coordinate  value  is  constant  and  equal  to 
the  negative  of  the  focal  point  C-f)  of  the  camera.   Using 
these  parameters  the  equation  becomes 


x-xo 

y-yO 
-f 


=    K 


a    b    c 
d    e    f 

g   h    i 


XA-XL 
YA-YL 

ZA-ZL 


C3.1) 


where  K  is  a  scale  factor.   From  this  transformation  the 
following  equations  are  obtained 

x-xo  -  K  [a  CXA-XLD  +  b  C YA-YL )  +  c  CZA-ZLDD   C3.2D 
y-yo  =  K  [d  CXA-XL3  +  e  CYA-YL)  +  f  C ZA-ZL 3 :   C3.3D 
-f   ■  K  [g  CXA-XL)  +  h  CYA-YL)  +  i  CZA-ZL^3   C  3  .  4  ) 
Dividing  the  last  equation  into   the  first  two  and  rear- 
ranging yields 
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xo  -  F 


a  CXA-XL)  +  b  CYA-YL)  +  c  CZA-ZL) 


g  CXA-XL)  +  h  CYA-YL.}  +  i  CZA-ZL) 


C3.5D 


yo 


d  CXA-XL}  +  5  CYA-YLD  +  c  CZA-ZL) 


g  CXA-XLD  +  h  CYA-YL)  +  l  CZA-ZL) 


C3.6) 


where  x  and  y  are  the  2D  image  plane  coordinates . CReF .  3, 
pp.  141-1423 

The  original  parameters  of  the  M  matrix  were 
calculated  from  the  w ,  0,  and  k  that  were  given  by  the  DMA 
with  the  original  photograph  .   These  represent  the  physical 
position  of  the  image  plane  in  relation  to  the  georect- 
angular  coordinate  system  at  the  time  the  picture  was  taken 
The  Focal  point  was  also  supplied  and  is  particular  to  the 
camera  that  was  used  to  take  the  original  photograph.   When 
the  synthesized  image  is  generated,  the  image  plane  is 
oriented  to  a  position  For  the  desired  viewing  angle,  which 
means  that  the  parameters  oF  the  M  matrix  will  change  and 
must  be  recalculated .   The  steps  used  to  determine  the 
desired  orientation  oF  the  image  plane  coordinates  will  be 
discussed  later  in  this  chapter. 

2 .   Reference  Image  to  Screen  Coordinate  TransFormation 
Qnce  the  x  and  y  image  coordinates  have  been 
calculated  they  are  translated  to  I  and  J  values  oF  the 
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original  image.   This  is  accomplished  using  the  affme 
transform  which  represents  a  2D  into  ED  coordinate  transfor- 
mation.  The  equation  that  accomplishes  this  is  derived  from 


C  3  .  7  ) 
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where  CI  and  CE    are  the  values  that  translate  the  image 
plane  origin  to  the  I  and  J  coordinate  system  origin.   They 
are  calculated  by  setting  I  and  J  to  0  and  salving  for  x  and 
y  .   To  get  the  desired  transformation  of  the  image  coor- 
dinates to  I  and  J  coordinates  the  inverse  transform  is 
taken.  CRef.3,  p. 5333 
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Again  the  original  J,  k,  1,  m,  CI,  and  CE    were 
supplied  by  the  DMA  with  the  original  image.   Equation  3.B 
represents  any  general  ED  into  ED  transformation,  and  it  is 
used  both  to  translate  the  original  image  plane  coordinates 
into  the  I  and  J  coordinates  of  the  reference  image  and  to 
transform  the  image  plane  coordinates  of  the  synthesized 
view  into  screen  coordinates  of  IA,  and  JA  . 

The  screen  coordinates  were  assigned  the  parameters 
IA,  which  represents  the  columns,  and  JA,  which  represents 
the  rows.   The  point  IA.JACl.l)  is  mapped  to  the  upper  left 
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corner  of  the  screen  and  IA, JAC51E, 51H3  is  the  lower  right 
corner  . 

The  screen  coordinates  represent  the  location  of  a 
pixel  within  a  frame.   To  convert  I  and  J  original  coor- 
dinates to  IA  and  JA  screen  coordinates  one  needs  to  know 
the  particular  frame  one  is  working  in  and  compensate  for 
the  difference  in  the  starting  location  of  the  frame 
coordinates  and  the  original  image  coordinates.   This  is 
accomplished  by  using  the  following   equations, 

IA  -  CI  -  I _ Frame)  C3.9) 

JA  -  C4339  -  J  .Frame  -  J)     C3.10D 

The  I  Frame  and  J  Frame  values  must  therefore  be 
given  to  determine  the  screen  coordinates  within  a  desired 
frame.   To  allow  flexibility  in  determining  which  frame 
image  would  be  used  to  extract   the  reference  image,  an 
interactive  input  of  the  IFrame  and  JFrame  coordinates  was 
appropriate . 

3 .   Sunthesized  Image  Plane  Rotation 

For  the  synthesized  views  that  were  generated,  the 
affine  transform  parameters  CI,  CE ,  J,  k,  1,  and  m   that 
would  map  the  newly  rotated  image  plane  into  the  screen 
coordinate  system  were  selected.   Since  the  image  planes  in 
the  synthesized  views  were  oriented  for  an  observer  looking 
north,  the  image  plane  coordinates  were  generated  by  cal- 
culating the  z-axis,  which  points  south,  from  two  georec- 
tangular  coordinate  vectors  calculated  from  two  different 
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terrain  data  points  along  the  same  longitude  line.   By 
taking  the  difference  between  the  X,  Y,  and  Z  coordinates  a 
third  vector  was  formed  that  defined  the  image  plane  z-axis. 
The  image  plane  y-axis  was  calculated  by  using  only  one  of 
the  terrain  data  points  used  to  calculate  the  z-axis.  By 
taking  the  negative  of  the  georectangular  coordinates  of  the 
terrain  data  point,  a  vector  is  produced  that  points 
downward  through  the  center  of  the  earth.  This  was  the  image 
plane  y-axis.   Once  the  z  and  y  axes  are  calculated,  the 
cross  product  of  y  cross  z  was  used  to  calculate  the  x-axis. 
Figure  3.1  demonstrates  the  resulting  image  plane. 


th  Center 


Fig.  3.1   Synthesized  Image  Plane  Coordinates 

This  image  plane  coordinate  system,  from  the 
observers  perspective  at  CO,0,-f),  would  have  the  x-axis 
pointing  left,  the  y-axis  painting  down,  and  the  z-axis 
pointing  directly  toward  the  observer.   The  screen 
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coordinates  have  the  IA  axis  pointing  right  and  the  JA  axis 
pointing  down.  Therefore,  the  new  values  of  CI  and  C2  with 
IA,  and  JA  equal  to  0  were  as  follows: 

CI  =  Maximum  assigned  x-image  value    (3.11) 
C3  -  0  C3.1BD 

which  aligned  the  image  plane  x,y(o,a)  to  the  screen 
coordinates  of  IA.JACO.Cn.   The  J,  k,  1,  and  m  values  of  the 
transformation  matrix  were  selected  to  scale  the  image  plane 
to  the  screen. 

Having  determined  the  three  synthesized  image  plane 
coordinate  vectors  in  terms  of  georectangular  vectors,  it  is 
relatively  easy  to  generate  the  M  matrix  parameters  using 
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were  S,  R,  and  T  are  the  georectangular  coordinate  vectors 
of  the  x,  y,  and  z  axes  of  the  rotated  image  plane.   This 
defines  the  new  transformation  matrix  that  will  be  used  to 
generate  the  synthesized  view  by  mapping  the  georectangular 
vectors  of  the  terrain  data  into  the  new  image  plane. 
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B.   TERRAIN  DATA 

The  Digital  Terrain  Elevation  Data  (DTED)  supplied  by 
the  DMA  came  as  a  rectangular  grid  of  elevation  data  points. 
Each  elevation  data  point  was  recorded  as  an  integer  value 
in  meters  above  sea  level.   If  a  particular  elevation  point 
was  unknown,  it  was  assigned  a  value  of  -327E7  to  assure 
that  it  would  not  be  confused  with  any  valid  elevation  data 
points.   The  rectangular  terrain  grid  listed  elevation 
points  every  one  second  of  a  degree  change  in  latitude  and 
longitude.   The  southwest  corner  of  the  terrain  grid  was 
defined  as  being  located  at  37*  22'  47"  N.  latitude  and 
-122'  05'  03"  U .  longitude.   From  this  reference  point  the 
elevation  data  was  laid  out  in  210  rows  by  233  columns.   The 
rows  represented  lines  of  constant  latitude  and  the  columns 
lines  of  constant  longitude.   With  this  information  the 
northeast  corner  of  the  terrain  grid  was  calculated  as  being 
located  at  37'  26'  17"  N.  latitude  and  -122  01' 04"  W. 
longitude . 

1 .   Data  Uerification 

The  first  problem  was  to  compare  how  accurately  the 
elevation  data  matched  up  with  the  original  image  data. 
This  required  taking  specific  elevation  points  that  were 
known  to  match  specific  image  paints,  then  translating  the 
georectangler  coordinates  of  those  elevation  points  to  the  I 
and  J  coordinates  for  comparison  to  the  original  image.   The 
w,  0!,  and  k  for  the  M  matrix  to  transform  georectangular  to 
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image  plane  coordinates  and  the  parameters  For  the  affine 

transform  from  image  plane  to  original  image  I  and  J 
coordinates  were  supplied  by  the  DMA  with  the  original  image 
data  . 

To  select  specific  elevation  points  for  comparison 
required  finding  a  method  of  distinguishing  unique  elevation 
patterns  that  could  match  specific  objects  in  the  image. 
The  technique  used  was  to  visually  display  the  elevation 
data  as  an  image.   The  elevation  image  file  was  produced  by 
assigning  grey  scale  values  to  each  elevation  data  point 
with  the  lower  elevations  receiving  the  darker  shades  and 
the  higher  elevations  the  lighter  shades.   When  the  eleva- 
tion image  was  displayed  as  shown  in  Figure  3.E,  a  distinct 
highway  pattern  emerged  from  which  three  intersections  could 
reasonably  be  distinguished.   The  three  intersection  points 
Cshown  as  a,  b,  and  c  in  Fig.  3.E)  correspond  to  elevated 
roads  that  crossed  one  another  in  the  original  image. 

Dnce  the  elevation  points  were  selected  for  com- 
parison, the  approximate  row  and  column  of  the  elevation 
data  corresponding  to  the  center  of  the  intersections  was 
determined.   From  the  reference  point  of  the  terrain  grid 
there  is  a  1  second  change  in  latitude  and  longitude  for 
each  row  and  column  which  allowed  the  calculation  of  the 
latitude  and  longitude  for  each  of  the  three  reference 
elevation  points.   The  latitude  and  longitude  far  each  point 
was  converted  to  georectangular  coordinates  using  a 
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conversion  program,  then  transformed  to  I  and  J  coordinates 
using  the  provided  DMA  parameters  as  explained  earlier. 
Each  of  the  three  elevation  points  mapped  to  within  10 
pixels  of  the  original  image  in  both  the  I  and  J  coordinate 
directions.   This  equates  to  less  than  1  second  error  in 
latitude  and  longitude  which  was  deemed  precise  enough  to 
establish  the  correlation  between  the  terrain  and  image 
data  . 


Fig.  3.2   Elevation  Image 

5 .   Elevation  Line  Drawing 

By  selecting  a  smaller  area  of  the  DTED  data,  a  more 
distinct  picture  could  be  studied.   An  intersection  used  to 
verify  image  and  elevation  correlation  was  selected  to  be 
used  as  the  reference  image.   A  smaller  rectangular  set  of 
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terrain  data  that  would  map  to  the  intersection,  plus  a 
small  section  of  the  surrounding  area,  was  extracted  From 
the  reference  terrain  grid.   This  smaller  set  of  terrain 
data  was  taken  from  rows  71  through  73  and  columns  172 
through  183  of  the  terrain  model . 

To  verify  that  this  set  of  elevation  points  would 
appear  like  the  desired  reference  image,  a  line  drawing 
illustrated  in  Figure  3.3  was  created  using  a  commercial 
graphics  program  called  MDUIE.BYU  CRef.  4).   This  program 
can  generate  a  connected  line  drawing  from  a  set  of  eleva- 
tion data  and  allows  rotation  as  well  as  magnification  of 
the  drawing.   To  assure  a  sharp  visual  contrast,  the 
elevation  data  was  magnified  by  a  factor  of  10  and  the 
drawing  rotated  to  a  useful  viewing  angle. 

The  results  were  not  as  desired  which  gave  an  early 
indication  that  the  resolution  of  the  terrain  data  may  not 
be  adequate  enough  to  generate  a  synthesized  view  that  would 
be  a  close  approximation  of  the  reference  image.   This  was 
further  verified  when  the  synthesized  view  was  produced  at  a 
later  time.   An  artificial  set  of  image  and  elevation  data 
was  used  to  verify  that  the  transformation  program  func- 
tioned as  desired  before  generating  synthesized  views  of  the 
original  reference  image.   The  artificial  elevation  points 
mapped  into  the  artificial  image  exactly  way  as  the  original 
elevation  data  would  map  to  the  original  reference  image. 
The  artificial  reference  image,  shown  in  Figure  3.4, 
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Fie.  3.3   Elevation  Line  Drawing 


.  3.4   Artificial  Reference  Image 


was  of  a  small  square  structure  resting  on  tap  of  a  much 
larger  square  structure.   Selecting  a  large  abject  Far  an 
artificial  reference  image  would  minimize  the  effects  of  the 
law  resolution  of  the  elevation  data.   This   produced  more 
reasonable  synthesized  images  that  demonstrated  the  perspec- 
tive transformation  more  accurately.   The  results  of  the 
transformation  will  be  discussed  further  on  in  this  chapter 
after  considering  some  of  the  algorithms  used  to  generate 
the  synthesized  view. 

C.   SPECIFIC  ALGORITHMS 

1 .   Image  Referencing  Algorithm 

Due  to  the  resolution  mismatch  between  the  reference 
image  and  the  terrain  data,  a  smaller  subset  of  the  terrain 
model  was  selected.   This  would  allow  the  transformation  of 
smaller  and  more  distinct  images  that  could  show  the  effects 
of  the  program  better  .   The  desired  terrain  elevation  points 
are  extracted  from  the  larger  reference  terrain  grid  by 
interactively  selecting  the  proper  rows  and  columns  that 
define  those  particular  points.   The  program  allows  up  to  a 
50  by  50  array  of  elevation  points  to  be  extracted.   This 
size  was  chosen  to  limit  the  time  required  to  generate  a 
synthesized  view.   From  the  rows  and  columns,  the  latitude 
and  longitude  is  tabulated  for  each  point,  and  is  then 
converted  to  georectangular  coordinates.   The  georectangular 
coordinates  are  written  to  a  file  in  row  order  from  left  to 
right  and  from  bottom  to  top.   The  first  set  of  coordinates 
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match  the  southwest  elevation  point,  and  the  last  set  of 
coordinates  the  northeast  elevation  point. 

The  georectangular  coordinates  of  each  terrain  data 
point  is  translated  to  I ,  and  J  original  image  coordinates 
using  the  camera  parameters  supplied  by  the  DMA .   They  are 
further  processed  to  IA  and  JA  screen  coordinates  using  the 
I_ Frame  and  J  Frame  of  the  desired  reference  image.   The  IA 
and  JA  coordinates  derived  from  the  elevation  points  will 
nou  correspond  to  the  IA  and  JA  coordinates  of  the  reference 
image.   Any  four  adjacent  elevation  points  will  define  a 
rectangle  which  is  divided  into  two  triangular  panels  by 
defining  a  line  connecting  two  opposite  corners.   Dnce  the 
triangular  panels  of  the  elevation  points  are  defined  in 
terms  of  IA  and  JA  coordinates,  the  corresponding  pixel  grey 
scale  values  of  the  reference  image  that  fall  within  the 
same  screen  coordinates  of  the  triangular  panel  are  col- 
lected and  averaged.   As  seen  in  the  example  of  Figure  3.5, 
the  calculated  average  grey  scale  value  is  placed  into  a 
file  along  with  the  three  specific  elevation  points  that 
make  up  the  triangular  panel .   This  procedure  is  repeated 
until  the  entire  terrain  grid  has  been  processed.    Each 
triangular  panel  represents  a  sample  or  average  grey  scale 
value  of  the  reference  image. 

When  the  latitude,  longitude,  and  elevation  of  the 
new  observation  point  is  input  into  the  program,  a  new  XL, 
YL,  and  ZL  is  calculated  that  corresponds  to  the  new 
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perspective  center.   As  explained  earlier,  the  image  plane 
is  rotated  to  the  desired  viewing  angle,  which  changes  the  M 
matrix  parameters  as  well. 
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Fig.  3.5   Triangular  Panel  Data  File  Structure 

Uith  these  new  parameters  for  the  perspective  transformation 
equations,  the  georectangular  coordinates  of  the  terrain 
grid  are  once  again  run  through  the  perspective  transfor- 
mation, and  the  new  IA  and  JA  coordinates  are  calculated  for 
each  point.   These  new  I A  and  JA  coordinates  represent  the 
transformed  elevation  points   as  seen  from  the  new  observer 
location.   The  next  step  is  to  take  the  same  three  elevation 
points  that  farmed  a  triangular  panel  in  the  original 
transformation,  map  them  into  the  synthesized  image  file 
using  the  new  IA  and  JA  coordinates,  and  then  fill  them  in 
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with  the  assigned  grey  scale  value.   The  mapping  and  Filling 
process  is  explained  in  the  next  section. 

2 .   Polygon  Fill  and  Hidden  Surface  Elimination 

In  the  perspective  transformation  of  the  terrain 
grid  into  the  synthesized  image  plane,  some  oF  the  tri- 
angular panels  became  partially  or  entirely  hidden  by  other 
panels.   To  determine  which  pixels  will  be  visible  and 
therefore  written  to  the   synthesized  image  and  which  pixels 
are  hidden,  the  z-buFFer  algorithm  was  used. 

Uhen  the  new  observer  location  is  input  into  the 
program,  the  XL,  YL ,  and  ZL  georectangular  coordinates  are 
tabulated.   Using  the  georectangular  coordinates  previously 
calculated  For  each  oF  the  terrain  grid  elevation  points, 
the  distance  or  depth  From  the  observer  location  to  the 
elevation  data  points  can  be  calculated  by  the  Following 
equation 

Depth  -  Square  root C CXL-X) 2  +  CYL-Y32  +C2L-2)2J    C3.14D 

The  depth  oF  the  pixels  within  a  triangular  panel 
were  calculated  by  using  the  normalized  plane  equation  that 
defines  the  plane  oF  the  triangular  panel.   The  normalized 
plane  equation  in  3D  space  is  given  by 

ax  +  by  +  cz  -  -1        C3.15) 
To  solve  For  the  coefficients  a,  b,  and  c,  the  three  eleva- 
tion points  that  specify  a  triangular  panel  are  used,  and 
the  x,  y,  and  z  are  replaced  with  IA,  JA ,  and  the  Depth  of 
each  elevation  point.   If  the  three  points  are  CIA1,  JA1 , 
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Depthl),  CIAS,  JAS,  DepthE:,  and  CIAS,  JA3,  Depth3),  then  in 
matrix  form  we  have  the  following 

IA1  JA1  Depthl 
IA3  JAS  DepthS 
IA3    JA3    Depth3 

Solving  for  the  coefficients  a,  b,  and  c  uie  have  CRef  .  3,  p. 
E0B3 


a 

-1 

b 

= 

-1 

C3.16) 

c 

-1 



.... 

a 

b 
c 


IA1  JA1 
IAS  JAS 
IA3     JA3 


Depthl 
DepthE 
Depth3 


-1 


-1 
-1 
-1 


C3.17) 


The  inverse  of  the  elevation  point  matrix  is 
determined  by  calculating  the  adjoint,  which  is  the  trans- 
pose of  the  cofactor  matrix,  and  multiplying  each  term  by 
the  reciprocal  of  the  determinant  CRef.  5,  p.  A-153 .   The 
depth  of  each  pixel  within  a  tranformed  triangular  panel  can 
now  be  determined  by  using  the  IA  and  JA  of  the  pixel  in  the 
following  equation. 

Depth  -  -  CI  +  a  CIA)  +  b  CJA))/c      C3.1BD 
Each  triangular  panel  will  have  different  a,  b,  and  c 
coefficients,  therefore  the  pixels  within  each  triangular 
panel  will  have  a  different  depth  than  those  from  another 
panel  . 

The  z-buffer  is  a  two  dimensional  array  that  is  the 
same  size  as  the  synthesized  image  of  51E  by  51E  pixels. 
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When  the  IA  and  JA  coordinates  of  a  pixel  is  determined,  the 
depth  is  calculated  and  then  compared  to  any  previously 
written  depth  in  the  z-buffer  at  that  same  coordinate 
location.   IF  the  depth  is  smaller,  it  means  that  the  pixel 
is  closer  to  the  observer  and  would  cover  the  previously 
written  pixel.   The  depth  oF  the  new  pixel  will  then  replace 
that  oF  the  previous  pixel.  The  assigned  grey  scale  value, 
which  is  determined  From  the  projected  triangular  panel  in 
which  the  new  pixel  is  contained,  is  written  to  the  syn- 
thesized image  at  the  calculated  IA  and  JA  coordinate 
location.   IF  the  new  pixel  depth  is  larger  then  the 
previous  pixel  depth,  that  means  the  new  pixel  is  Farther 
away  From  the  observer  and  would  be  covered  by  the  previous 
pixel.   In  this  case  no  action  is  taken  and  the  next  pixel 
is  processed.   This  procedure  is  continued  until  all  the 
pixels  within  each  translated  triangular  panel  has  been 
processed.   CReF.  2,  pp.  265-2673 

The  IA  and  JA  coordinates  oF  the  pixels  in  the 
synthesized  image  are  calculated  From  each  transformed 
triangular  panel  using  an  active  edge  list  (AEL),  J  Bucket, 
and  Frame  buFFer .   The  IA  and  JA  screen  coordinates  oF  the 
three  elevation  points  that  deFine  a  triangular  panel  are 
used  to  generate  the  parameters  oF  the  AEL .   Three  lines  are 
Formed  that  connect  each  oF  the  three  translated  elevation 
points  and  are  identiFied  as  line  1,  2,  and  3.   A  line  is 
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defined  between  any  tuia  paints  from  which  you  can  determine 
the  IA  associated  with  the  maximum  JA  of  the  two  points  by 
IA_INCPT  -  C  IA  value  of  the  maximum  JA  )        C3.19) 
How  much  IA  changes  for  a  one  step  change  in  JA  is  given  by 


DELTA  IA  = 


IACPomt  S3  -  IACPoint  1) 


JACPoint  3.1  -  JACPoint  13 


C3.E0) 


which  is  the  inverse  slope  of  the  line  between  the  two 
paints.   The  span  of  JA  between  two  paints  is  given  by 

DELTA__JA  =  JACPoint2)  -  JACPoint  1)  C3.21D 

For  each  line,  the  IA  coordinate  that  corresponds  to 
the  maximum  JA  value  of  the  line,  the  amount  IA  changes  for 
each  one  unit  step  of  JA,  and  the  total  span  of  JA  is 
determined  and  stared  into  the  AEL .   After  the  line  para- 
meters are  placed  into  the  AEL,  the  line  identification 
number  is  put  into  the  J  Bucket ,  which  is  a  one  dimensional 
array  of  size  512,  at  the  same  location  as  the  maximum  JA  of 
the  line.   In  this  way  the  J_Bucket  acts  as  a  pointer  to  the 
lines  stored  in  the  AEL.   If  a  previous  line  number  has 
already  been  written  to  that  location,  then  the  line  ident- 
ification number  of  the  new  line  is  placed  into  the  AEL  and 
is  linked  to  that  line  already  residing  in  the  J_.Bucket .   An 
example  is  shown  in  Figure  3.5  where  line  #1  is  referenced 
by  the  J  Bucket  and  line  #2  is  linked  to  line  #1.   The  IA 
and  JA  coordinates  of  the  pixel  locations  to  be  mapped  to 
the  synthesized  image  are  contained  within  the  three  lines 
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whose  parameters  are  contained  in  the  AEL .   If  a  line  is 
horizontal,  the  IA  and  JA  coordinates  of  that  line  are 
located  between  the  other  two  lines  and  therefore  does  not 
need  to  be  written  to  the  AEL .  CRef.  2,  pp. 75-783 


J  Bucket 


AEL  For  Line  1 


Line 
#1 


AEL  For  Line  2 
Fig.  3.6   Active  Edge  List  Storage 

The  J_Bucket  is  scanned  from  its  maximum  to  minimum 
coordinate  value.   If  the  J  Bucket  contains  a  line  pointer, 
the  parameters  for  that  line  are  retrieved  from  the  AEL  as 
well  as  those  of  any  line  it  is  linked  to.   Since  we  have 
defined  a  triangle  there  will  always  be  two  lines  to  work 
with.   From  the  parameters  of  the  two  lines,  the  IA  coor- 
dinates that  fall  between  them  is  calculated  for  each  JA 
scan  line.   A  scan  line  is  all  the  columns  CIA  coordinates) 
along  a  particular  row  CJA  coordinate).   As  the  JA  coor- 
dinate is  decremented  by  one,  from  maximum  to  minimum,  the 
J Bucket  is  checked  to  see  if  a  new  line  has  been  added,  and 
then  the  corresponding  IA  for  each  line  is  determined  for 
that  JA  value.   Those  IA  values  and  any  that  fall  between 
them  are  mapped  to  the  frame  buffer . 
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The  frame  buffer  is  a  512  by  512  array  that  is 
initialized  to  0 .   As  the  IA  and  JA  values  are  mapped  into 
the  frame  buffer  the  coordinate  locations  matching  the  IA 
and  JA  values  are  changed  from  0  to  1 .   This  process 
continues,  and  each  time  the  JA  scan  line  is  decremented, 
the  DELTA  JA  parameters  for  each  of  the  two  lines  are  also 
decremented  and  the  J Bucket  is  checked.   If  the  J  ^Bucket 
contains  another  line  pointer,  then  the  line  it  references 
will  replace  the  line  whose  DELTA  JA  has  decreased  to  0. 
When  finished,  the  frame  buffer  will  have  recorded  the  IA 
and  JA  coordinates  for  every  pixel  location  within  the 
translated  triangular  panel  .   The  frame  buffer  is  then 
scanned,  and  if  a  particular  IA  and  JA  coordinate  location 
contains  a  value  of  1,  then  the  depth  for  a  pixel  located  at 
those  same  coordinates  is  calculated  and  compared  to 
previously  tabulated  depths  in  the  z-buffer  as  explained 
earlier . 

This  process  is  known  as  a  polygon  fill  routine  that 
utilizes  the  z-buffer  algorithm  for  hidden  surface  elimina- 
tion.  If  any  part  of  a  triangular  panel,  after  going 
through  the  perspective  transformation,  should  map  outside 
the  synthesized  image  coordinate  boundary,  the  entire  panel 
is  discarded  and  the  next  panel  is  processed.   This  is  not  a 
satisfactory  solution  and  could  have  been  corrected  by 
implementing  a  clipping  routine  that  would  allow  partial 
triangular  panels  to  be  mapped  to  the  synthesized  image. 
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However,  due  to  time  constraints,  a  clipping  routine  was  not 
incorporated  into  the  3D  transformation  program. 

D.   SUnriARY 

Using  the  artificial  reference  image  from  Figure  3.4  and 
an  artificial  terrain  grid  that  mapped  to  the  same  image, 
two  synthesized  views  were  generated.   The  first  synthesized 
view  shown  in  Figure  3.7  was  from  an  observation  point 
located  at  37°  22'  30"  N.  latitude,  -122°  01*  59"  W. 
longitude,  and  an  elevation  of  110  meters.   Figure  3.8 
exhibits  the  next  synthesized  view  that  was  produced  from  an 
observation  point  of  37*  20'  0"  N.  latitude  while  main- 
taining the  same  longitude  and  observer  height  as  before. 
This  simulates  viewing  the  object  from  further  away .   As 
expected  of  a  perspective  view  the  object  appears  smaller. 
The  near  view  also  demonstrates  the  perspective  relationship 
by  the  apparent  taper  from  front  to  back.   A  third  perspec- 
tive view  is  depicted  in  Figure  3.9  from  close  in  and  at  a 
higher  elevation.   The  observer  location  was  from  37"  23' 
30''  N.  Latitude,  -122*  01'  59''  Ul .  Longitude,  and  a  height 
of  150  meters.   This  demonstrates  the  visual  effect  of  not 
displaying  partial  triangular  panels  in  the  image  and  why  a 
clipping  routine  is  needed. 

The  actual  reference  image  is  shown  in  Figure  3.10  from 
which  the  synthesized  view  displayed  in  Figure  3.11  was 
generated.   In  comparing  the  reference  image  to  the  syn- 
thesized view  there  is  little  resemblance.   A  closer 


synthesized  view  is  seen  in  Figure  3. IS  which  gives  a 
partial  representation  of  the  reference  image.   The  lack  of 
detail  in  both  shading  and  the  depiction  of  features  can  be 
attributed  to  the  poor  resolution  of  the  terrain  data  Cap- 
proximately  25  meter  resolution)  as  compared  to  that  of  the 
reference  image  data  CI  meter  resolution) . 

To  improve  the  quality  of  the  synthesized  view  the 
terrain  data  must  be  of  a  higher  resolution.   Having  more 
elevation  data  points  over  a  given  area,  the  fewer  number  of 
reference  image  pixels  one  must  collect  and  average  for  a 
triangular  panel.   The  synthesized  image  will  then  maintain 
a  closer  approximation  to  the  various  shades  of  the  refer- 
ence image.   The  physical  shape  of  an  object  also  suffers 
from  poor  terrain  data  resolution.   The  perspective  trans- 
formation forms  straight  lines  between  translated  elevation 
points.   This  causes  object  distortion  if  the  elevation 
points  do  not  fall  exactly  along  the  boundaries  of  the 
object.   As  an  example,  if  a  square  box  in  the  reference 
image  had  only  one  elevation  point  defined  on  its  surface  at 
the  center  of  the  box,  then  the  synthesized  image  can  not 
reproduce  the  corners  and  edges  of  that  box,  and  it  would 
appear  distorted.   For  these  reasons  the  closer  the  resolu- 
tion of  the  terrain  data  to  the  resolution  of  the  reference 
image,  the  closer  the  synthesized  view  resembles  the  refer- 
ence image  in  shading  and  shape. 
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Fig.  3.7   Artificial  Synthesized  Image  1 


Fig.  3.8  Artificial  Synthesized  Image  2 
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Fig.  3.3  Artificial  Synthesized  Image  3 


Fig.  3.10   Reference  Image 
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Fig  3.11  Synthesized  Image  1 


Fig.  3.12   Synthesized  Image  E 
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The  3D  transformation  program  was  developed  to  allow 
tracing  the  Flow  of  data  easily.   Much  of  the  data  was 
written  to  files  so  that  it  could  be  printed  out  and 
studied.   This  method  of  data  storage  required  certain  files 
to  be  read  many  times  as  the  synthesized  image  was  gener- 
ated.  The  program  execution  would  have  been  faster  if  the 
data  had  been  stored  in  arrays  that  are  easily  passed 
between  the  various  modules.   Another  consideration  that 
would  increase  speed  would  be  to  decrease  or  eliminate  the 
interactive  input  required.   This  could  be  accomplished  by 
extracting  the  desired  terrain  data  into  a  file  and  sepa- 
rating the  associated  reference  image  before  program 
execution.  This  would  require  longer  set  up  time  but  would 
decrease  the  amount  of  data  manipulation  required  by  the 
program  and  increase  its  overall  speed. 
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iu.   CONCLUSIONS 

A .   GENERAL 

The  3D  computer  image  transformation  from  a  photographic 
image  was  a  difficult  task  to  achieve  satisfactorily.   The 
main  objective  was  to  develop  a  program  that  takes  a  grey 
scale  photographic  image,  a  set  of  elevation  data  points 
defined  over  that  image,  and  generates  a  rotated  synthesized 
perspective  view.   This  goal  was  realized,  but  some  areas 
still  need  improvement.   The  quality  of  the  synthesized  view 
is  Judged  on  how  well  the  grey  scale  values  of  pixels  match 
those  of  the  original  image  and  how  closely  the  translated 
objects  resemble  the  desired  structure. 

The  resolution  of  the  terrain  model  as  compared  to  the 
resolution  of  the  reference  image  data,  extensively  affects 
the  synthesized  image  quality.   Depending  on  the  contents  of 
the  reference  image,  the  required  resolution  of  the  terrain 
model  will  vary.   If  the  image  is  of  open  country,  the 
distance  between  elevation  points  may  be  large  and  the 
result  may  not  suffer  unacceptable  degradation  in  the  syn- 
thesized view.   If  the  image  is  of  a  city  that  has  many 
small  distinct  objects  such  as  buildings,  then  the  resolu- 
tion of  the  elevation  points  must  be  higher.   Given  a  set  of 
elevation  data  points,  the  question  to  be  resolved  is  how  to 
make  the  synthesized  pixels  relate  to  the  reference  image 
better,  and  thus  improve  the  quality  of  the  synthesized 
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views7   This  question  and  alternative  ways  to  improve  speed 
and  program  flexibility  will  be  discussed. 

B.   GREY  SCALE  CORRELATION 

There  are  a  few  possible  methods  to  improve  the  shading 
of  the  synthesized  view  to  better  match  that  of  the  original 
image.   The  translated  triangular  panels  form  very  distinc- 
tive lines  or  boundaries  between  areas  of  different  shades. 
It  may  be  desirable  to  blend  these  boundaries  by  sampling 
the  pixels  along  both  sides,  replacing  those  pixels  with  an 
averaged  value.   Another  possibility  would  be  to  implement  a 
Gouraud  or  Phong  shading  algoritm  CRef.  3,  pp.  323-3303 . 
This  would  help  smooth  the  shade  transition  across  the 
boundary  and  result  in  a  more  smooth  appearance. 

Another  method  to  improve  grey  scale  correlation  would 
be  tc  develop  a  way  to  divide  the  triangular  panels  into 
smaller  triangular  areas  before  referencing  the  original 
image.   By  having  smaller  triangular  panels,  smaller  areas 
of  the  reference  image  are  sampled  resulting  in  a  better 
approximation  in  the  synthesized  view. 

Only  the  left  image  of  a  stereo  photographic  pair  of 
images  was  used  in  this  study  .   The  right  image  may  contain 
grey  scale  information  of  surfaces  hidden  in  the  left  image 
view  or  vice  versa.   By  sampling  both  left  and  right  images 
and  complementing  them,  some  ambiguities  may  be  resolved. 
These  suggestions  will  increase  the  number  of  calculations 
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required  to  generate  the  synthesized  view  but  may  improve 
the  quality  of  the  synthesized  view. 

C.   PROGRAM  SPEED  AND  FLEXIBILITY 

Although  speed  was  not  a  prime  consideration  in  this 
study,  it  would  be  desirable  to  generate  synthesized  views 
as  quickly  as  passible.   To  determine  which  subroutines 
consumed  the  mast  CPU  time  the  program  was  monitored  while 
running  the  artificial  reference  data  set.   Table  1  shows 
the  results  of  the  CPU  time  used  and  the  percentage  of  the 
total  time  far  each  subroutine  called  by  the  main  program. 
If  a  subroutine  calls  another  subroutine  that  time  is 
included  in  the  calling  routines  CPU  time.   The  subroutines 
that  required  interactive  input  were  not  measured.   There- 
fore, the  total  CPU  time  cannot  be  measured  precisely. 
However,  the  results  obtained  represent  reasonable  estimates 
and  demonstrate  where  the  focus  should  be  on  improving 
overall  speed. 

Some  suggestions  have  already  been  discussed  on  how  to 
improve  the  speed  of  the  program,  but  ather  methods  also 
exist.   Preconditioning  the  input  data  to  elemmate  interac- 
tive input  and  using  arrays  instead  of  files  were  two 
methods  already  presented.   The  z-buffer  is  accessed  several 
times  during  synthesized  view  generation.   If  the  z-buffer 
could  be  implemented  in  hardware,  the  time  required  for 
processing  of  polygon  fill  and  hidden  surface  elimination 
routines  would  be  improved. 
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TABLE 

1 

CPU 

TIME  CONSUMPTION 

CIN  SECONDS) 

SUBROUTINE 

CPU  TIME 

PERCENTAGE  C5sD 

TER  CROP 

14.  11 

9.051 

READIMAGE 

3.  23 

2.072 

TER  INTRP 

0.04 

0.0B5 

REAL  EL 

0.07 

0.045 

TER._DMS 

O.BS 

0.443 

in  REFIJ 

0.33 

0.E50 

IM__REFAUG 

1  .10 

0.70B 

AFFIN 

0.04 

0.026 

NEW  I J 

O.BO 

0.513 

NODE  ..DEPTH 

0.B3 

0.404 

FILL 

134. BO 

B6.4B7 

TOTAL 

155.3 

100.0 

The  frame  buffer  is  used  and  accessed  in  two  different 
areas  of  the  program  for  each  triangular  panel  translated. 
It  may  be  passible  to  eliminate  this  buffer  by  processing 
each  pixel  as  its  IA  and  JA  coordinates  are  determined, 
instead  of  staring  that  information  into  the  frame  buffer 
for  later  processing.   Other  polygon  fill  routines  such  as 
seed  fill  algorithms  or  using  fence  registers  may  be  faster 
than  the  edge  fill  routine  used  in  this  program  CRef  .  2,  pp. 
BO-BBJ  . 

For  program  flexibility  it  uould  be  desirable  to  be  able 
to  adjust  the  image  plane  of  the  synthesized  view  to  any 
desired  angle.   The  program  limits  the  image  plane  to  a 
northern  direction.   To  make  the  synthesized  image  plane 
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adjustable  would  require  developing  a  method  for  rotating 
the  image  plane  coordinates  in  terms  of  the  georectangular 
coordinates.   This  mould  improve  the  3D  transformation 
program  and  extend  its  usefulness.   Another  program  improve- 
ment would  be  the  incorporation  of  a  clipping  routine  to 
improve  the  appearance  of  partial  synthesized  views  as 
discussed  in  previous  sections. 

The  ideas  used  to  develop  this  program  were  contrived 
from  fundamental  concepts.   Many  areas  were  discussed  that 
could  improve  or  enhance  the  basic  implementation  of  the 
program.   The  results  that  were  obtained  are  encouraging  and 
could  easily  be  used  as  the  basis  for  further  study. 


APPENDIX  A 

PR  El  2  PA'"1  sunnARY 


I  .   PROGRAM  TRAN  3  E 

a.   F  u  n  c tiers  eerfrrmed 

Takes  an  image  and  elsvatian  data  File  ard  extracts 
the  des.rsd  regie"-  far  transformation.   Accepts  interactive 
input  cf  a  new  cbserver  location  and  generates  a  s^nt^es.zsd 
zieuj  da  making  calls  tc  ./arieus  subroutine  noddies.   This  is 
the  Ta:n  cart  of  the  program. 


e  . 


Passes  parameters  bet-ee^  subroutines. 

^n  iMPGES  file  that  contains  the  synthesized  vie_. 

Galling  routines 

\cne  . 

Called  routines 

INPUT:  Cbtams  the  interactive  input  ef 

the  elevation  a^d  mags  file  name: 


and  desired  reus  and  ooi 


!Sl    _^!  =  — 


-i   p-.'f  r; 


•'  'a' 


"Ep  EFG1-1 


:-    Cata  points 
*~  ~ — n  *~  ^  e  r3^a,-ap^a  sle-atio- i  c r  — 
E  :tracts  the  dssiret  elevation 
points  cf  the  area  to  be  trans- 
lated . 


cp 


READIMAGE 


TEE  INTRP 


PE£L  EL 


TER  CHS: 


Reads  the  reference  image  into  ar 

arraL,  called  in^GE. 

This  rQLitins  collects  and  av ara^ss 


<-hi 


ins  extracted  elevation  poi 


ar.d 


P  E  F  I  J  : 


assigns  that  value  to  a^y  unknown 


elevation  data  points. 
Unites  the  extracted  elevation 
points  into  a  file  called  ZFIL.D^T 
as  real  values  in  meters  ate -ye  saa 
level . 

Determines  the  latitude  and  lon- 
gitude of  the  elevation  points  then 
converts  them  to  georectangular 
coordinates  and  stores  them  in  a 
file  called  XYZ.DAT. 
Converts  the  georeotangular  rccr- 


ir    EEEPUE 


AFFI 


and  J  A  screen  coordinates. 

Constructs  the  file  NODE. DAT  that 
contains  the  three  elevation  points 
and  grey  scale  '/alee  that  make  up  a 
triangular  canel . 

Determines  the  affme  transform 
coefficients  utilized  in  the 
transform  of  the  synthesized  image 


CQ 


plans    cccrciirate! 
lates  . 


:□  scree: 


i—.  i »  -  -: 


CBS  LDC:  Accepts  ths  interactive  input  cf 

the  latitude,  longitude,  and  height 

in  meters  above  sea  level  cf  the 

desired  observer  location. 
NEW  I  J:  Computes  ths  neuj  I  A  and  J  A  scrsen 

coordinates  cf  the  elevation  points 

For  the  synthesized  view. 
NODE  CPTH:        Calculates  the  distance  from  each 

s  1  e v aticn  point  tc  the  observer 

location . 
FILL:  Determines  the  hidden  Surfaces  and 

generates  the  synthesized  image  in 

the  file  iriAGES.DAT. 
Pontine  parameters 
IENDN:  The  number  cf  reus  cf  the  extracted 

elevation  data  points. 
IENDM:  The  number  cf  columns  cf  the 

extracted  elevation  data  points. 


a.       SUBROUTINE  INPUT 

a,   Functions  performed. 

Accepts  the  interactive  input 
_mage  data  files,  as  well  as  the  infer 
extract  the  tssirsd  elevation  points - 


cf  ths  elevation  a^d 

mation  needed  tc 


b.   Input  (interactive) 


IROU: 


LRTUJ: 


I  Cui_ 


E'  C  T  '  c  . 

I M  F I L  E  : 
I  FRAME : 

JFRAHE : 


Minimum  rou  of  the  elevation  area 

desired . 

Ma:-:  imum  rcuj  cf  the  elevation  area 

desired . 

Minimum  column  value  cf  the  desired 

elevation  area. 

Maximum  column  value  cf  the  desired 

elevation  area. 

The  reference  elevation  file  name. 

The  reference  image  file  name. 

The  I  Frame  value  cf  the  reference 

image . 

The  J  Frame  value  cf  the  reference 


image . 
c    Dutput 

Same  as  the  interactive  input 

d.  Tailing  routines 
IRAN  3  D  . 

e.  Tailed  routines 
None  . 

f  .   Routine  parameters 
None  . 

3.   SUBROUTINE  TER  TRTP 

a.   Functions  Performed 


SI 


smal ler 
IELEUS  . 
b  . 


Reads  the  original  terrain  grid  and  constructs  a 

grid  of  the  desired  elevation  points  in  the  array 


Input 

LRGUJ: 
T  rnr  . 

LCGL: 
ELFIL: 

Output 

IELEUS:  An  array  containing  the  extracted 

elevation  points. 
Galling  routines 
TPPN  3  D  . 
Called  routines. 
None  . 

Routine  parameters 
IELEU1:  An  array  containing  the  entire 

elevation  data  set. 


Minimum  rcu  of  the  desired 

elevation  area . 

naximum  rcu  of  the  desired 

elevation  area. 

riinimum  column  value  of  the  desired 

elevation  area . 

Maximum  column  value  cf  the  desired 

elevation  area. 

The  file  containing  the  reference 

elevation  data . 


a  Xi 


Counters 


4.       SUBROUTINE    REAOiriAGE 

a.   Functions  performed 


IMAGE 


Reads  the  reference  image  into  an  array  called 


T  r-|  —  |  ■  f- 


IflFIL : 


The  file  containing  the  reference 
image  pixel  grey  scale  data. 


I"!  AGE: 


An  array  containing  the  pi:-: el  grey 
scale  values  of  the  reference 
image . 


Galling  routines 


T?AM  3  D . 
Galled  routines 


Routine  parameters 

Counters 


P  =mH   T 


IK  a no  IL: 


5.   SUBROUTINE  TER  INTRP 

a.  Functions  performed 

Determines  the  average  elevation  over  the  extracted 
elevation  area  and  assigns  that  elevation  to  any  unknown 
elevation  points. 

b.  Input 

IENDN:  Number  cf  extracted  elevation  reus 

Number  of  extracted  elevation 
columns . 


CT 


IELEU2: 

Gutput 
IELEU2 : 


ar  arraL,  cf  the  extracted  Elevation 
data  . 

ar:aL,  ccptai".°g  the  extracted 
elevation  data  :.=^y    unknown  data 
3  o  i  n  t  ^  nave  ^app  ace'  -mQ-^  *-  ^  a 

average  elevation  of  the  area.."' 
Galling  routines 
TRAN  3  G  . 

Galled  routines 

Mone  : 

Routine  parameters 

N  and  M :  Gcunters . 

NEXT:  Count  of  the  number  of  elevation 

points  processed. 
IEL:  Summation  cf  the  elevations. 

I  AUG:  The  average  elevation  of  the  area. 


SUBROUTINES  PEAL  EL 

a    ^" j~ct 1  crs  nsrf crTsd 

Greates  a  real  file  called  ZFIL. 


GAT    cf    the    e: 


•t-earl 


-l-l- 


CS  . 


b.        Input 
IENDN 

t  r r  irjM 


T[na    nufT,!zBr    cf    extracted    sls'aJ 

reus  . 

TKe  number  of  extracted  eleva' 

ccium^s  . 


E4 


IELEU2: 

uUtput 

ZFIL.DAT 


An  array  cf  the  extracted  elevation 
data  points. 

A  File  ccrta:":rg  the  real  v a  1  u  e  □  f 
the  extracted  elevation  data 
points . 


.a  i  i 


1  i  *— i  —r 

'a 


r3utir,SS 

T?AN  3  D  . 
Galled  routines 
None  . 

Routine  parameters 

PELEU:  A  real  array  cf  the  e: 

elevation  paints. 


traoted 


2MS 


a.  Functions  performed 
[Idvsrts  each  extracted  elevation  data  point  to  its 

equivalent.   It  uses  the  fact  that  eaoh  elevation  point 
represents  a  one  second  change  in  latitude  cr  longitude  from 
the  next  point,  starting  from  the  Southwest  corner  reference 
ele;  a:::n  located  at  3"°  32'  4""  N.  Latitude  and  -122°  05' 
03"  LJ  .  Longitude.   The  DMS  is  converted  to  X,  Y  and  2 
gecrectangu i ar  coordinates  and  stored  in  the  XYZ.DAT  file. 

b.  Input 

IRC'jJ:  The  minimum  row  cf  the  extracted 

elevation  data  points. 


S3 


LRQU 


ICOL: 


IENDH: 


!FIL  .DA' 


Output 
XYZ.EAT 


The  maximum  rou  oF  extracted  eleva- 

tian  data  paints. 

The  minimum  column  value  aF  the 
extracted  elevation  data  paints. 
The  maximum  column  value  of  the 
extracted  elevation  data  paints. 
The  total  number  of  reus  cf  the 
extracted  elevation  data  points. 
The  File  containing  the  real  values 
oF  the  extracted  elevation  data 
points . 

Q  File  containing  the  gecrec- 
tangular  coordinates  cF  the  ex- 
tracted elevation  data  paints. 


Calling  routines 


IRAN  3  D. 

Tailed  routines 
CHS2XYZ . 

Routine  parameters 
IL  and  JL: 


V ,  and  Z 


£>T' 


LATH,  and 


Counters 

Gecrsctangular  coordinates  oF  an 

sis',  at  ion  paint  . 

The  latitude  in  degrees,  minutes, 

a^d  seconds  or  an  elevation  point 


ES 


LDND,  LDNM, 

and  LQNS :         The  longitude   in  degrees,  minutes, 

and  seconds  of  an  elevation  point. 
HEIGHT:  The  real  elevation  in  meters  above 

sea  level  of  an  elevation  point. 
R  L  A  T  *'  and  P  I_  A  T  S  :  The  rsf srsncs  latitude  in  m m_i t a c 

and  seoonds  of  the  reference  eleva- 

t i □  n  point  located  at  r o u  1,  column 

1 

RLDNn  AND  RLQNS:  The  reference  longitude  in  minutes 

and  seconds  of  the  reference  eleva- 
tion pcmt  located  at  rou  1  ,  column 


S.   SUBROUTINES  in  REFIJ 

a.  Functions  performed 

Calculates  the  I  and  J  reference  image  coordinates 
and  converts  them  to  the  IA  and  JA  screen  coordinates  for 
each  extracted  elevation  data  point  and  puts  them  into 
arrays  IA,  and  JA  . 

b.  Input 

IFRAME:  The  I  Frame  value  of  the  reference 


JFRAME 


"Mnri  • 


image . 

t j-j pa    Frari6  vsiue  of  the  reference 

image . 

The  number  cf  reus  in  the  extracted 

elevation  data. 


IENDN 


XYZ  .  DA' 


Gutput 
IA: 


JA: 


The  number  of  cclumns  in  the  ex- 
tracted elevation  data. 


Tata  File  containing  the  g  s  c  r  e  c  - 


tangular  coordinates  of"  the 
extracted  elevation  data  paints. 

Pn  array  containing  the  I  A  screen 
coordinate  values  of  the  extracted 
elevation  data  points. 
An  array  containing  the  JA  screen 
coordinate  values  of  the   extracted 
elevation  data  points. 


Calling  routines 
TPAN  3  D  . 
Galled  routines 
PROJECT,  XY2IJ. 
Routine  parameters 


OMEGA: 


PH 


KAPPA: 


XQ  and  YG 


Rctaticn  cf  the  original  image 
plane  about  x-axis  in  radians. 
Rctaticn  cf  the  original  image 
plane  about  the  y-axis  in  radians 
Rctaticn  cf  the  original  image 
plane  about  the  z— axis  in  radians 
Cffset  cf  the  principal  point  Frcr 
the  I  PP. 


E3 


XI,  Yl,  and  21 

Al,  A2,  31,  B2 

CI,  and  CE: 


p  o  c  u  s 


xima 


Y  I  MA 


I  am 


The  gscrsctanguiar  cocrdinatas  c  f 
the  observer  location. 

The  given  affine  transfers  para- 
meters for  the  transform  of  the 
original  image  plane  coordinates  to 
I  and  J  original  image  screen 
coordinates . 

The  camera  focal  length. 
The  x-axis  value  of  the  elevation 
point  in  image  plane  coordinates. 
The  Y-axis  value  of  the  elevation 
point  in  image  plane  coordinates. 
The  original  reference  image 
screen  coordinates. 


SUBROUTINE  DMSEXYZ 

a.  Functions  performed 

Converts  the  CMS  latitude  and  longitude  data 
°nd  Z  georectangular  coordinates. 

b .  I n^ut 

LATD,  LATH,  and 

LQTS:  The  latitude  m  degrees,  min 

and  seconds  of  the  extracted 


utes 


.5  . 


I M     =  T-l  ■ 


LCNS 


HEIGHT 


The  lonoit—de  in  decrees   fninLt3s 
and  ssccnds  c  f  the  extracted 
elevation  data  points. 
The  height  in  meters  above  sea 


data  point 


■  i    *  i 


and  Z.: 


—  —  wLUj.'  a  _S±:=   _L 


TS-fpH   != 


] a 1 i i n g  routines 


TEE  DHS   GEE  LEE 


Tailed  routines 

None  . 

Routine  parameters 


PHI 


LAMDA 


N ,  E  .  SQUARE , 

and  A: 


P  3  ^  *  ~.  p '  . 


PI 


Pngle  in  radians  from  the  Z-axis  of 
the  georectanguiar  coordinate 

system  . 

Angle  in  radians  from  the  X-axis  of 

the  georectanguiar  coordinates 
system . 

Given  parameters  used  in  calcu- 
lating the  georectanguiar  cocr- 
dinatas  of  the  elevation  pci^is. 
Number  of  radians  per  degree. 
Number  of  radians  in  a 


CI: 

C2: 


number  of  degrees  in  a 
Number  of  minutes  in  a  degree 
Number  of  seconds  in  a  decree 


SUEPQUTINE  P  P  G  _T  E  C  T 


a  . 


mcturs  performed 


Converts  the  X,  Y ,  and  2  georectanguiar  cccrdi^a-as 
of  the  extracted  elevation  data  paints  ta  the  ::-imags  and 


y   i  .  1 1  a  y  = 


t . 


the  image  piai 


OMEGA :  x-axis  rotation  of  the  image  plane 

in  radians. 
PHI:  y-axis  rotation  of  the  image  plane 

in  radians . 
KAPPA:  z-axis  rotation  of  the  image  plane 

in  radians . 
XQ  and  YD:        The  offset  of  the  principal  point 

from  the  IPP. 
XI,  Y 1 ,  and  Zi:   The  georectanguiar  coordinates  of 

the  observer  location. 


i    4  i 


and  Z: 


FDCUS: 


The  georectanguiar  coordinates  of 


the  extracted  elevation  dat; 
pcints . 

The  camera  focal  length. 


xinfi 


YIMA: 


The  image  plane  x-axis  coordinate 
value  of  the  extracted  elevation 
data  points . 

The  image  plans  y-axis  coordinate 
value  of  the  extracted  elevation 


joints 


Tailing  routines 

in  pefi j . 

Tailed  routines 
None  . 

Routine  parameters 
Till,  HIE,  M 1 3 , 
MEi,  nEE,     ME3, 
m31,  M3E,  M33: 


The  parameters  of  the  M  matrix. 
The  denominator  of  the  transforma- 
tion equations. 


1 .  SUBROUTINE  XVEIJ 

a.   Functions  performed 

Converts  the  image  plane  coordinates  of  the  ex- 
racted  elevation  data  points  to  the  I  and  J  original  imagt 
ccrc^nates . 

h  .        Input 


X  I  HA  : 


The  image  plane  x-axis  coordinate 
values  of  the  extracted  elevation 
data  points . 


IE 


Yin  a 


fll,  AE,  El,  EE 
CI,   EE: 


Output 
I  and  J 


The  image  plane  y-axis  coordinate 
values  of  the  extracted  elevation 
data  points . 

The  aff ins  transform  parameters  to 
map  the  image  plane  coordinates  to 
the  reference  image  I  and  J 
coordinates . 

The  reference  image  coordinate 
values  for  the  extracted  elevation 
data  points . 
d.   Calling  routines 
in  REFIJ,  NEW  .  I  J  . 
s.   Called  routines 

None  . 
f  .   Routine  parameters 

DENQM:  The  denominator  cf  the  aff ire 

transform . 

IE.  SUBROUTINES  in-REFAUS 

a.   Functions  performed 

Constructs  the  file  NODE . DAT  that  contains  the 
elevation  points  and  reference  grey  scale  value  that  make  up 
a  triangular  panel. 


Inpu' 


I A  and  J  A 


The  arrays  that  contain  the  re- 
ference image  screen  coordinates  of 
the  extracted  elevation  data 
pci n ts . 

The  array  that  contains  the 
reference  image  grey  scale  values. 
The  number  of  rows  in  the  extracted 
elevation  data. 
The  number  of  columns  in  the 
extracted  elevation  data. 


IM^GE: 
IENDM: 
IENDN: 

Output 

NODE. DAT:         The  file  that  contains  the  eleva- 
tion paints  that  make  up  a  tri- 
angular panel  and  its  associated 
reference  grey  scale  value. 

Calling  routines 

TRflN  3  D  . 

Called  routines 

None  . 

Routine  parameters 

SLOPE:  The  slope  of  the  diagonal  line  that 

separates  the  rectangle,  defined  by 
four  elevation  data  pcints,  into 
two  triangular  panels. 

VINT:  The  Y  intercept  of  the  diagonal 

line . 


IY,  II,  and  N 

I  GREY  1: 


I  GREYS: 


L,  LI,  and  IP: 
NODE  A,  NODE  5 
NODE  C,  and 
NODE  D: 


I  COUNT 1 : 


ITOTi 


I  COUNTS 


Counters . 

The  averaged  reference  grey  scale 

value  of  the  triangular  panel  below 

the  diagonal  line. 

The  averaged  referenced  grey  scale 

value  cf  the  triangular  panel  above 

the  diagonal  line. 

Counters . 


The  numerical  designation  cf  the 

four  elevation  points  that  make  up 

the  rectangle  that  is  divided  into 

two  triangular  panels. 

The  count  of  the  number  of  pixels 

averaged  in  the  triangular  panel 

belou  the  diagonal  line. 

The  summation  of  the  pixel  grey 

scale  values  averaged  in  the 

triangular  panel  belou  the  diagonal 

1  me  . 

The  ccunt  of  the  number  of  pixels 

averaged  in  the  triangular  panel 

above  the  diagonal  line. 

The  summation  of  the  pixel  grey 

scale  values  averaged  in  the 


-*c; 


triangular  panel  above  the  diagonal 
line . 


13.  SUBROUTINE  EXT  QPIEN 

a.  F unctions  performed 

Determines  the  fl  matrix  parameters  used  m  the 
rotation  of  the  image  plane  for  the  synthesized  view.   This 
program  defines  the  three  image  plane  coordinates  in  terms 
of  georeotanguiar  coordinates. 

b.  Input 

IENDM:  The  number  of  rous  in  the  extracted 

elevation  data  points. 

IENDN  The  number  of  columns  in  the 

extracted  elevation  data  points. 


Output 

(111,  il  1 2  ,  M 1 3  : 

MSI,  M22,  M23: 

M31,  f13S,  M33: 

Calling  routines 
New  _I  J  . 

Tailed  routines 

N'cne  . 

Routine  parameters 


First  row  coefficients  of  the 

transformation  matrix. 

Second  row  coefficients  of  the 

transformation  matrix. 

Third  row  coefficients  of  the 

transformation  matrix. 


MAGN  X: 


HAGN 


:  iHurj  ^  : 


(.  V,  and  Z: 


i    *  i 


X  CGRD: 


V  CGRD: 


:grd 


X  UECX,  X  .UECY, 

and  X  UECZ: 


Y  UECX,  Y  UECY 


and  Y  UECZ 


The  magnitude  of  the  synthesized 
view  image  plane  x-axis. 
The  magnitude  cf  the  synthesized 
vieu  image  plane  y-axis. 
The  magnitude  cf  the  synthesized 
view  image  plane  z-axis. 
Georectangular  coordinates  cf  the 
extracted  elevation  data  points. 
Array  that  stores  the  georec- 
tangular  X  coordinates  of  the 
extracted  elevation  data  points. 
Array  that  stares  the  georsc- 
tangular  Y  coordinates  of  the 
extracted  elevation  points. 
Array  that  stores  the  georec- 
tangular Z  coordinates  cf  the 
extracted  elevation  data  points. 

The  X,  Y,  and  Z  georectangular 
coordinates  of  the  synthesized 
image  plane  x-axis. 

The  X,  Y,  and  Z  georectangular 
coordinates  of  the  synthesized 
image  plane  y-axis. 
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2 UECX,  Z  UECY 

and  Z  UECZ : 


IP  a"d  IP 


The  X,  Y,  and  Z  georsctangular 

-Dardinatss  af  the  synthesized 

image  plane  z-a x i s . 

The  tctal  number  of  extracted 

elevation  points. 

Counters . 


14.  SUBROUTINE  AFFIN 

a.  Functions  performed 

Assigns  or  calculates  the  coefficients  to  be 
utilized  in  the  affme  transform  from  image  coordinates  to 
screen  coordinates  of  the  synthesized  view. 

b.  Input 
None  . 

c .  Output 

<=rl,  AS,  El,  BE 

and  CI,  C2 :       The  affme  transform  coefficients 

for  the  synthesized  view  image 
plane  coordinates  to  screen 
coordinates . 

d.  Calling  routines 
IRAN  3  D . 

e.  Called  routines 
None 

f .  Routine  parameters 


XII1A  HAX 


yiha  riAx 


J  MAX 


Assigned  maximum  image  plane  :•: 

coordinate  . 

Assigned  maximum  image  plane  y 

coordinate . 

Assigned  maximum  IA  screen 

coordinate . 

Assigned  maximum  JA  screen 

coordinate  . 


15.  SUBROUTINE  GES  LGC 

a.   Functions  performed 

Calculates  the  neu  observer  location  georectangulai 
coordinates  from  the  interactive  input  of  the  desired 
latitude,  and  longitude,  and  height.   This  routine  also 
assigns  the  focal  length  for  the  synthesized  view, 
b  ■   Input  '"interactive;' 

LATH  and  LATS:    The  minutes  and  seconds  of  the 

latitude  of  the  neu  observer 
location . 
LQNM  and  LDNS :    The  minutes  and  seconds  of  the 

longitude  of  the  neu  observer 


HEIGHT: 


location 


'e  altitude  in  meters  above  sea 


level  for  the  neu  observer  loca- 


icn  . 


d  . 


XI,  Yl ,  and  21:   The  gscrectangular  coordinates  of 

neu  observer  location. 

FOCUS:  Assigned  Focal  length  For  the 

synthesized  view. 

Calling  routines 

TRAM  3  D  . 

Called  routines 

DMS2XYZ . 


F .   Routine  parameters 
None 

IE.  SUBROUTINE  NEU  IJ 

a.   Functions  perFormed 

Calculates  the  neu  IA,  and  JA  screen  coordinates  oF 
the  extracted  elevation  data  points  using  the  neu  observer 
location  data . 
b  .   Input 

XI,  Yl  and  21:    Gecrectangular  coordinates  oF  the 

neu  observer  location. 
FOCUS:  Assigned  Focal  length  For  the 

synthesized  view. 
Al,  AE,  El,  EE, 

and  CI,  CR-.  The  aFFine  transform  ocsFFicients  . 

I ENDfl :  The  number  oF  reus  oF  the  extracted 

elevation  data  points. 


BO 


IENDN: 

Output 
IA: 


JA: 


The  number  of  columns  of  the 
extracted  elevation  data  points. 

The  array  containing  the  elevation 
data  IA  screen  coordinates  cf  the 
synthesized  view. 

The  array  containing  the  elevation 
data  JA  screen  coordinates  cf  the 
synthesized  view. 


f  . 


Calling  routines 
TRAN  3  D  . 
Called  routines 
XY2IJ,  EXT  OR  I  EN. 
Routine  parameters 
nil,  ni2,  and 

ni3: 


M21,  M22,  and 

H23: 

1131,  H32,  and 

M33: 

XD,  and  YO: 

XI  HA: 


First  row  coefficients  of  the 
transformation  matrix. 

Second  rou  coefficients  of  the 
transformation  matrix. 

Third  row  coefficients  cf  the 

transformation  matrix. 

The  offset  of  the  principal  point 

to  the  I  PP. 

Image  plane  x  coordinate. 
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YiriA: 
DENDM 

ITOT: 
IR: 


Image  plane  y  coordinate. 

Denominator  of  the  transformation 

matrix  . 

Total  number  of  extracted  elevation 

data  points  . 

Counter  . 


17.  SUBROUTINE  NDDE  DPTH 

a.  Functions  performed 

Calculates  the  distance  from  each  transformed 
elevation  data  point  to  the  neu  observer  location. 

b.  Input 

XI,  Yl,  and  21:   The  georectangular  coordinates  of 

the  neu  observer  location. 

IENOM:  The  number  of  raus  in  the  extracted 

elevation  data  points. 

IENDN:  The  number  of  columns  in  the 

extracted  elevation  data  points. 


Output 
DEPTH: 


Array  of  the  distances  from  the 
transformed  elevation  data  points 
to  the  neu  observer  location. 


Calling  routines 

TRAN  3  D  . 

Called  routines 

None  . 

Routine  parameters 
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IR: 

I  T  G  T  : 


Gaunter . 

Tata!  number  of  extracted  si eve' 

data  points . 


19.  SUEHOUT !*•.'£  FILL 

a.   Functions  performed 

Determines  the  hidden  surfaces  using  a  z-buffer 
algorithm.   The  depth  is  compared  for  each  pixel  at  a 
specified  screen  coordinate  tc  determine  if  it  is  written  to 
the  synthesized  image. 
t.   Input 

IP:  Prray  containing  the  IP  screen 

coordinates  of  the  transformed 
elevation  data  pcmts. 
JA :  Array  containing  the  JA  screen 

coordinates  of  the  transformed 
elevation  data  points. 
CEPTH:  Array  of  the  distance  from  the 

transformed  elevation  data  points 
tc  the  neu  observer  location. 
IMAGE:  Array  used  to  ccnstruct  the 

synthesized  image. 
IENDH:  The  number  of  rows  of  the  extracted 

elevation  data  points. 
IENDN:  The  number  cf  columns  of  the 

extracted  elevation  data  points. 
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XYZ . DAT : 


c  . 


Output 
IMAGES: 


File  of  the  georectangular  coor- 
dinates of  the  elevation  data 
points . 

A  file  containing  the  synthesized 
image . 


Calling  routines 

TPAN  3  C  . 

Called  routines 

FRAME  FIL. 

Routine  parameters 

Z  EUFF:  The  z-buFFer 

XI   Yi   Zl   XE 

Y'2   ZE   v3   Y3 

and  Z3:  The  screen  coordinates  and  depth  cF 

the  three  elevation  points  that 

deF_ne  a  triangular  panel. 
Z  EPTH:  The  distance  cF  a  pixel  within  a 

triangular  panel  to  the  new 

observer  location. 
Cll,  C1E,  C13, 
C21,  C22,  C23, 
C31,  C32,  C33:    The  coFactcrs  cF  the  plane  equation 

matrix . 
GET:  The  determinate  af  the  plane 

equaticn  matrix. 
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A  COEF,  B  COEF, 
and  C  CDEF: 

I  GREY: 

NDDE  A,  NODE  E, 
NODE  C: 


I  M  I  N  : 
I  MAX: 
J  J1IN: 

J  MAX: 

FRAME: 
IR,  I,  J, 
K,  and  L: 

I  PLANES: 


The  coef f icients  of  the  plane 
equation . 

Grey  scale  value  of  a  triangular 
panel . 

Numerical  dssignaticn  cf  the  three 

elevation  data  paints  that  make  a 

triangular  panel. 

Minimum  IA  screen  coordinate  of 

the  three  elevation  points. 

Maximum  IA  screen  coordinate  cf  the 

three  elevation  points. 

Minimum  JA  screen  coordinate  cf  the 

three  elevation  points. 

Maximum  JA  screen  coordinate  cf 

the  three  elevation  points. 

The  frame  buffer  . 

Counters . 

Total  number  of  triangular  panels 
constructed  from  the  extracted 
elevation  data  points. 


es 


19.  SUBROUTINE  FRAME  FIL 

a  .   Functions  performed 

Constructs  an  edge  list  from  the  three  transformed 
gi Hvs^iori  do i n  t s .   t]-|oco  edoes  form  the  hcundar',as  ^ r~ rm    a 
polygon  fill  routine.   The  pi :< els  that  are  determined  to 
fall  uiithin  the  transformed  triangle  are  marked  in  the  frame 
buffer  . 

b .   I nput 

NODE  A,  "iGOE  E, 

and  NODE  C:       The  numerical  designation  of  the 

three  elevation  data  paints  that 
make  a  triangular  panel  . 
Array  that  contains  the  IA  screen 
coordinates  of  the  transformed 
elevation  data  points. 
Array  that  contains  the  JA  screen 
coordinates  of  the  transformed 
elevation  data  points. 
Maximum  JA  screen  coordinate  of  the 
three  elevation  data  paints. 
Minimum  JA  screen  coordinate  of  the 
three  elevation  data  paints. 

The  frame  buffer  array  . 


IA 


JA 


J  MAX  : 

J  MIN: 

Output 

FRAME: 

Calling  routine 

FILL. 


8E 


e.  Called  routine 
None  . 

f .  Routine  parameters 

I  INCPT:  The  matching  IA  coordinate  cf  the 

maximum  JA  paint  of  a  line. 

DELTA  I:  The  amount  Ifi  changes  Far  a  one 

step  change  in  JA . 

DELTA  J:  The  total  JA  span  cf  a  line. 

AEL :  Array  containing  the  parameters  cf 

the  three  lines  cf  transformed 
triangular  panel . 

XNDDE1  and 

XNODE  2:  Designates  the  number  of  I A  coor 

dmates  between  two  lines  for  a  given 

JA. 

DX:  Indicator  used  to  determine  the 

direction  in  which  the  IA  coor- 
dinates are  counted . 

NODE:  Array  containing  the  numerical 

designation  of  the  three  elevation 
data  points . 

NDDE1,  and 

NQDE2 :  Designators  for  determining  the 

elevation  point  that  contains  the 
highest  JA  v alue  between  two 
paints . 


N  HIGH  and 
N  LDU: 


IT,  I U ,  IP 
and  IS: 
J  5UCKET: 

II  and  12: 


ICNT  and  LCNT 


Used  with  NDDE1  and  NDDE2  Far 

determining  the  highest  JA  value 
betuieen  tuo  paints. 

Counters . 

Array  that  contains  the  line  number 

designators  referencing  the  PEL. 

Used  to  determine  the  number  of  IA 

coordinates  to  be  written  to  the 

frame  buffer . 

Counters  far  the  J  BUCKET. 
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APPENDIX  B 

3D- TRANS FORMAT I ON  PROGRAM  LISTING 

PROGRAM  TRAN_3_D 
C 

C        THIS  PROGRAM  TAKES  AN  IMAGE  AND  ELEVATION  FILE  AND 
C        CONSTRUCTSTHE  REFERENCE  IMAGE  AND  ELEVATION  FILE. 
C        FROM  THESE  FILES  A  SYNTHESIZED  IMAGE  IS  PRODUCED. 
C 

CHARACTER  ELFILE*13 , IMFILE*13 

BYTE  IMAGE(  512,512) 

INTEGER  IELEV2( 50,50) , IA( 2500) ,JA( 2500) 

I NTEGER  I ROW , LROW , I COL , LCOL , I ENDN , I ENDM , 

INTEGER  I FRAME ,  JFRAME 

REAL  XI , Yl ,Z1, FOCUS , DEPTH(  2500) 

REAL  Al , A2 , Bl , B2 , CI , C2 , OMEGA, PHI , KAPPA 
C 

CALL  INPUT( IROW,LROW, I COL, LCOL, ELFILE, IMFILE, 
I FRAME, JFRAME) 

OPEN(UNIT=l,FILE=ELFILE,STATUS=,OLD" ) 

OPEN(UNIT=2,FILE='ZFIL.DAT' ,STATUS='NEW' , 

ACCESS='DIRECT' , RECORDSIZE=128, MAXREC=512 ) 

OPEN(UNIT=3,FILE='XYZ.DAT' , STATUS=' NEW' , 

ACCESS=' SEQUENTIAL'  , FORM= *  FORMATTED '  ) 

OPEN( UNIT=4,FILE=IMFILE, STATUS=*OLD" , ACCESS=' DIRECT' , 
RECORDS I ZE=12 8 , MAXREC=5 12 ) 

OPEN( UN I T=2  0 , F I LE= ' NODE . DAT " , STATUS=  *  NEW '  , 

ACCESS=' DIRECT" , RECORDS  I ZE=128, MAXREC=5 12 ) 

OPEN(UNIT=21,FILE=' IMAGES.DAT' , STATUS= *  NEW* , 

ACCESS=' DIRECT' , RECORDS I ZE=128, MAXREC=5 12 ) 

CALL  TER_CROP( I ROW, LROW, I COL, LCOL, IELEV2 ) 

I ENDN=LROW- I ROW+ 1 

IENDM=LCOL-ICOL+l 

CALL  RE AD I MAGE (  IMAGE) 

CALL  TER_INTRP( IENDN, IENDM, IELEV2 ) 

CALL  REAL_EL( IENDN, IENDM, IELEV2 ) 

CALL  TER_DMS( I ROW, LROW, ICOL, LCOL, IENDM) 

CALL  IM_REFIJ( IA, JA, I FRAME, JFRAME, IENDM, IENDN) 

CALL  IM_REFAVG(  IA, JA, IMAGE, IENDM, IENDN) 

CALL  AFFIN(A1,A2,B1,B2,C1,C2) 

CALL  OBS_LOC( XI ,Y1,Z1, FOCUS) 

CALL  NEW_IJ( XI , Yl , Zl , FOCUS , Al , A2 , Bl , B2 , CI , C2 , 
IENDM, IENDN, IA, JA) 

CALL  NODE_DPTH( XI , Yl , Zl , DEPTH, IENDM, IENDN) 

CALL  FILL( IA, JA, DEPTH, IMAGE, IENDM, IENDN) 

CLOSE( 1) 

CLOSE(2) 

CLOSE( 3) 

CLOSE( 4) 

CL0SE(20) 
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CL0SE(21) 
END 
********************************************************** 

SUBROUTINE  INPUT( IROW,LROW, ICOL, LCOL, ELFILE, IMFILE, 

I FRAME, JFRAME) 

INITIAL  ROW  OF  DESIRED  AREA 
LAST  ROW  OF  DESIRED  AREA 
INITIAL  COLUMN  OF  DESIRED  AREA 
LAST  COLUMN  OF  DESIRED  AREA 

ELEVATION  DATA  FILE  NAME 

IMAGE  DATA  FILE  NAME 

I  FRAME  NUMBER 

J  FRAME  NUMBER 
CHARACTER  ELFILE*13 , IMFILE*13 
INTEGER  I ROW, LROW, ICOL, LCOL, IFRAME, JFRAME 


c 

c 

I  ROW  : 

c 

LROW  : 

c 

ICOL  : 

c 

LCOL  : 

c 

ELFILE 

c 

IMFILE 

c 

I FRAME 

c 

JFRAME 

35 


45 


55 


WRITE( 6,* 
WRITE(  6,* 
READ( 5,35 
WRITE( 6,* 
READ( 5,35 
WRITE( 6, * 
READ( 5,3  5 
WRITE( 6,* 
READ( 5,35 
FORMAT ( 13 
WRITE( 6,* 
READ(  5,45 
WRITE( 6,* 
READ( 5 , 45 
FORMAT( Al 
WRITE(  6,* 
READ( 5,55 
WRITE( 6,* 
READ(  5,55 
FORMAT(  14 
WRITE( 6,* 
RETURN 
END 


' INPUT 
' ENTER 
I  ROW 
' ENTER 
LROW 
1  ENTER 
ICOL 
1  ENTER 
LCOL 

' INPUT 
ELFILE 
1  INPUT 

IMFILE 


3) 


1  INPUT 
IFRAME 
*  INPUT 
JFRAME 


ELEVATION  AREA  DESIRED' 
MINIMUM  ROW  NUMBER  :  ' 

MAXIMUM  ROW  NUMBER  :  ' 

MINIMUM  COLUMN  NUMBER  :  ' 

MAXIMUM  COLUMN  NUMBER  :  ' 

THE  ELEVATION  DATA  FILE  NAME  :  ' 
THE  IMAGE  DATA  FILE  NAME  :  * 

THE  I  FRAME  NUMBER  OF  IMAGE  :  ' 
THE  J  FRAME  NUMBER  OF  IMAGE  :  ' 


*****  WAIT  APPROX.  1  MINUTE  FOR  SETUP****** 


C 

c 

c 
c 
c 
c 
c 


********************************************************** 
SUBROUTINE  TER_CROP( I ROW, LROW, ICOL, LCOL, IELEV2 ) 

THIS  SUBROUTINE  READS  THE  ORIGINAL  TERRAIN  GRID 
AND  CONSTRUCTS  A  SMALLER  GRID  OF  THE  ELEVATION 
POINTS  DESIRED. 

CHARACTER  SWLAT*8, SWLON*8, DELLAT*4, DELL0N*4, 
CHARACTER  COLS*4, R0WS*4 

INTEGER  IELEV1( 210,239) , IELEV2( 50,50) , IROW,LROW 
INTEGER  ICOL, LCOL, JV( 239) 

THE  VALUES  OF  JV, SWLON, SWLAT  DELLON, DELLAT, COLS 
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C      AND  ROWS  ARE  IGNORED. 

READ( 1,5) SWLON , SWLAT , DELLON , DELLAT , COLS , ROWS 
5      FORMAT( 1X,2(A8,2X),4(A4,2X) ) 

DO  20  M=l,238 
READ(  1,10)JV(M),(  IELEV1(N,M),N=1,210) 

io    format(  I6,2x,20i6/(8x,20i6) ) 

20     CONTINUE 

DO  40  N=IROW,LROW 
DO  30  M=ICOL/LCOL 
IELEV2(M+l-ICOL,N+l-IROW)=IELEVl(N,M) 
30      CONTINUE 
40     CONTINUE 
RETURN 
END 
C 

Q      ********************************************************** 

SUBROUTINE  READ I MAGE ( IMAGE ) 
C 

C        THIS  SUBROUTINE  READS  THE  IMAGE  DATA  INTO  AN  ARRAY. 
C 

BYTE  IMAGE(512,512) 
C 

DO  10  IR=1,512 
READ(  4/REC=IR)(  IMAGE(  IC, IR) , 10=1,512) 
10     CONTINUE 

RETURN 

END 
C 

Q      ********************************************************** 

SUBROUTINE  TER_INTRP( IENDN, IENDM, IELEV2 ) 
C 

C        THIS  SUBROUTINE  DETERMINES  THE  AVERAGE  VALUE  OF  THE 
C        ELEVATION  DATA  AND  ASSIGNS  THAT  VALUE  TO  ANY  UNKNOWN 
C        POINTS 
C 

INTEGER  IENDN, IENDM, IELEV2( 50,50) 
DATA  NEXT, I EL, IAVG/0,0,0/ 
C 

DO  20  N=l, IENDN 
DO  10  M=l, IENDM 
IF(  IELEV2(M,N).  EQ.  -32767)THEN 

GOTO  10 
ELSE 
NEXT=NEXT+1 
IEL=IEL+IELEV2(M,N) 
END  IF 
10      CONTINUE 
20     CONTINUE 

IAVG=IEL/NEXT 
C       CHANGE  UNKNOWN  ELEVATION  VALUES  TO  THE 
C        CALCULATED  AVERAGE. 
DO  40  N=l, IENDN 
DO  30  M=l, IENDM 
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IF( IELEV2(M,N). EQ. -32767)THEN 

IELEV2(M/N)=IAVG 
END  IF 
30      CONTINUE 
40     CONTINUE 
RETURN 
END 
C 

Q  ********************************************************* 

SUBROUTINE  REAL_EL( IENDN, IENDM, IELEV2) 
C 

C      THIS  SUBROUTINE  CREATS  A  REAL  FILE  OF  THE  ELEVATIONS 
C 

INTEGER  IELEV2( 50,50), IENDN, IENDM 
REAL  AELEV(239) 
C 

DO  20  N=l, IENDN 
K=0 

DO  10  M=l, IENDM 
K=K+1 

AELEV(K)=IELEV2(M/N) 
10      CONTINUE 

WRITE(2/REC=N)(AELEV(  INT), INT=1, IENDM) 
20     CONTINUE 
RETURN 
END 
C 

Q      ********************************************************** 

SUBROUTINE  TER_DMS( IROW, LROW, I COL, LCOL, IENDM) 
C 

C        THIS  SUBROUTINE  CONVERTS  EACH  ELEVATION  POINT  TO  ITS 
C        DMS  EQUIVALENT.  IT  USES  THE  FACT  THAT  EACH  ELEVATION 
C        POINT  REPRESENTS  A  ONE  SECOND  CHANGE  IN  LATITUDE 
C        OR  LONGITUDE  FROM  THE  NEXT  POINT. 

C        REFERENCE  POINT:  LAT.  037  DEG.  :  22  MIN.  :  47  SEC.  NORTH 
C  LON.  -122  DEG.  : 05  MIN.  : 03  SEC.  WEST 

C        OUR  ELEVATION  POINTS  STAY  WITHIN  THE  BOUNDS  OF  037 
C        DEGREES  NORTH  AND  -122  DEGREES  WEST,  SO  THESE  VALUES 
C        WILL  BE  ASSIGNED. 
C 

IMPLICIT  DOUBLE  PRECISION  (A-Z) 

INTEGER  IROW, LROW, I COL, LCOL, IENDM, IL, JL 

REAL  X , Y , Z , LATD , LATM , LATS , AELEV(  239) 

REAL  LOND,LONM,LONS, HEIGHT 

PARAMETER(RLATM=22.  ,RLATS=47.  ,RLONM=5.  ,RLONS=3.  ) 
C 

LATD=3  7. 

LOND=-122. 

DO  -20  I  L=  I  ROW,  LROW 
K=IL/60 
LATM=RLATM+K 
LATS=RLATS+( IL-K*60) 
IF(LATS.  GE.  60.  0)THEN 
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LATS=LATS-60. 0 
LATM=LATM+1. 0 
END  IF 
Il=IL-( IROW-1) 

READ(2/REC=I1)(AELEV(  INT), INT=1, IENDM) 
1=0 
DO  10  JL=ICOL,LCOL 
K=JL/60 
L0NM=RL0NM-K 
LONS=RLONS-( JL-K*60) 
IF(LONS.  LT.  0.  0)THEN 
LONM=LONM-l.  0 
LONS=LONS+60. 0 
END  IF 
LONM=-LONM 
L0NS=-L0NS 
1  =  1  +  1 

HEIGHT=AELEV(  I) 
CALL  DMS2XYZ( LATD, LATM, LATS, LOND, LONM, LONS, 

HEIGHT, X,Y,Z) 
WRITE( 3,99)X,Y,Z 
99        FORMAT(  3(  1X,F17.  7) ) 
10       CONTINUE 
20      CONTINUE 

ENDFILE(  3) 
REWIND( 3) 
RETURN 
END 
C 

SUBROUTINE  IM_REFIJ(  IA,  JA, IFRAME, JFRAME, IENDM, IENDN) 
C 

C        THIS  SUBROUTINE  CONSTRUCTS  THE  I  AND  J  COORDINATE 
C        DATA  FOR  EACH  ELEVATION  POINT  AND  THEN  CONVERTS  THEM 
C        TO  SCREEN  COORDINATES  AND  STORES  THEM  IN  ARRAYS 
C        I A  AND  JA. 
C 

IMPLICIT  DOUBLE  PRECISION  (A-Z) 

REAL  OMEGA , PHI , KAPPA , X , Y , Z , XI , Yl , Zl , XO , YO 

REAL  XIMA,YIMA  Al, A2 , Bl , B2 , CI , C2 , FOCUS 

INTEGER  I, J, IENDM, IENDN, IA( 2500 ), JA( 2500 ) 

INTEGER  IFRAME, JFRAME, ITOT 

DATA  OMEGA, PHI, KAPPA/.  8341764,-. 4563699,3. 0761254/ 

DATA  X0,Y0/0. 000002,0.0/ 

DATA  XI, Yl,Zl/-2693765.  9,-4304520.  4,3859018.  3/ 

DATA  Al,A2/20. 11323959,-6. 022849824/ 

DATA  Bl,B2/6. 016207940,20.  r0938801/ 

DATA  Cl,C2/-34954. 59484, -22566.  71593/ 

DATA  FOCUS/0.153197/ 
C 

I TOT= IENDN* IENDM 

DO  10  IR=1, ITOT 
READ( 3,5,END=20)X,Y,Z 
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5       FORMAT( 3(  1X,F17.  7) ) 

CALL  PROJECT( X, Y, 2 , OMEGA, PHI , KAPPA, XO, YO , XI , Yl , Zl , 

XIMA, YIMA, FOCUS) 
CALL  XY2IJ( XIMA, YIMA, I , J, AL,  A2,B1,B2,C1,C2 ) 

IA( IR)=I-IFRAME 
JA( IR)=4999-JFRAME-J 
10     CONTINUE 
20     REWIND(3) 
RETURN 
END 
C 

Q      ********************************************************** 

SUBROUTINE  DMS2XYZ( LATD, LATM, LATS, LOND, LONM, LONS, 

LONS, HEIGHT, X,Y,Z) 
C 

C        THIS  SUBROUTINE  CONVERTS  DMS  DATA  TO  X,Y,  AND  Z 
C        GEORECTANGULAR  COORDINATES. 
C 

IMPLICIT  DOUBLE  PRECISION  (A-Z) 

REAL  PHI , LAMDA,N,X, Y, Z, LATD, LATM, LATS 

REAL  LOND, LONM, LONS, HEIGHT 

PARAMETERS  I  =3.  14159265358793238) 

PARAMETER(C 1=180. , C2=60. ,C3=3600. ) 

PARAMETER( E_SQUARE=0.  006768658, A=6378206.  4) 
C 

RADIAN=PI/C1 

PHI=( LATD+LATM/C2+LATS/C3 ) *RADIAN 

LAMDA=( LOND+LONM/C2+LONS/C3 ) *RADIAN 

N=A/SQRT( 1-E_SQUARE*SIN(PHI)*SIN(PHI) ) 

X=(N+HEIGHT)*COS(PHI)*COS( LAMDA) 

Y=(N+HEIGHT)*COS(PHI)*SIN(LAMDA) 

Z=(N*( 1-E_SQUARE)+HEIGHT)*SIN(PHI) 

RETURN 

END 
C 

Q     ********************************************************* 

SUBROUT INE  PRO JECT( X , Y , Z , OMEGA , PHI , KAPPA , XO , YO , XI , 

Yl , Zl , XIMA, YIMA, FOCUS ) 
C 

C        THIS  SUBROUTINE  CONVERTS  THE  X,Y,Z  GEORECTANGULAR 
C        COORDINATES  TO  XIMA  AND  YIMA  WHICH  ARE  IMAGE 
C        PLANE  COORDINATES. 
C 

IMPLICIT  DOUBLE  PRECISION  (A-Z) 
REAL  Mll,M12,M13/M21,M22/M23,M31,M32,M33,DENOM 
REAL  XIMA, YIMA, X, Y, Z , OMEGA, PHI , KAPPA 
REAL  X0,Y0,X1,Y1,Z1, FOCUS 
C 

Mll=COS( PHI ) *C0S( KAPPA) 
M12=C0S( OMEGA) *SIN( KAPPA) + 

SIN( OMEGA) *SIN( PHI ) *C0S(  KAPPA) 
M13=SIN( OMEGA) *SIN( KAPPA) - 

COS( OMEGA) *SIN( PHI ) *COS(  KAPPA) 
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M21=-C0S(  PHI ) *SIN(  KAPPA) 

M22=COS(  OMEGA) *COS( KAPPA) - 

SIN( OMEGA) *SIN(  PHI ) *SIN( KAPPA) 

M23=S IN( OMEGA ) *COS(  KAPPA )  + 

COS(  OMEGA) *SIN(  PHI ) *SIN( KAPPA) 

M31=SIN(PHI) 

M32=-SIN(  OMEGA) *COS( PHI ) 

M33=COS(  OMEGA) *COS( PHI ) 
C 

DEN0M=M31*(X-X1)+M32*(  Y-Y1)+M33*( Z-Zl) 

XIMA=XO-FOCUS*(M11*(X-X1)+M12*( Y-Y1)+M13*( Z-Zl) )/DENOM 

YIMA=YO-FOCUS*(M21*(X-X1)+M22*(Y-Y1)+M23*(Z-Z1) )/DENOM 

XIMA=XIMA* 1000000 

YIMA=YIMA* 1000000 

RETURN 

END 
C 

Q      ********************************************************* 

SUBROUTINE  XY2IJ( XIMA, YIMA,  I , J, Al , A2 , Bl , B2 , CI , C2 ) 
C 

C        THIS  SUBROUTINE  TAKES  THE   IMAGE  POINTS  XIMA,YIMA 
C       AND  CONVERTS  THEM  TO  I , J  ORIGINAL  IMAGE  COORDINATES. 
C 

IMPLICIT  DOUBLE  PRECISION  (A-Z) 

REAL  XIMA, YIMA, Al , A2 , Bl , B2 , CI , C2 , DENOM 

INTEGER  I, J 
C 

DEN0M=A1*B2-B1*A2 

I=( ( XIMA-C1 ) *B2-( YIMA-C2 ) *B1 )/DENOM 

J=-( ( XIMA-C1 ) *A2-( YIMA-C2 ) *A1 ) /DENOM 

RETURN 

END 
C 

Q      ********************************************************* 

SUBROUTINE  IM_REFAVG( IA, JA, IMAGE, IENDM, IENDN) 
C 

C        THIS  SUBROUTINE  CONSTRUCTS  THE  FILE  THAT  IDENTIFIES 
C        THE  ELEVATION  POINTS  THAT  MAKE  UP  A  TRIANGULAR  PANEL 
C        AND  ITS  ASSOCIATED  SAMPLED  GREY  SCALE  VALUE. 
C 

IMPLICIT  DOUBLE  PRECISION  (A-Z) 
REAL  SLOPE ,YINT 
BYTE  IMAGE( 512,512) 

INTEGER  IA( 2500 ),JA( 2500), IY,M,N, IGREY1, IGREY2,L,L1 
I NTEGER  I ENDM , I ENDN , NODE_A , N0DE_B , N0DE_C , N0DE_D , I R 
INTEGER  IC0UNT1 , IC0UNT2 , IT0T1 , IT0T2 
C 

DO  90  IR=1, I ENDN- 1 
IL=( IR-1)* IENDM 
DO  80  N=1+IL,IENDM-1+IL 
NODE_A=N 
NODE_B=N+l 
N0DE_C=N+ IENDM 
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SLOPE=( JA(  NODE_C ) - JA( NODE_B ) ) * 1 . 0/ 

( IA(NODE_C)-IA(NODE_B) )*1.  0 
YINT=1.  0*JA(NODE_B) -SLOPE* IA(NODE_B) 
IT0T1=0 
ITOT2=0 
ICOUNT1=0 
ICOUNT2=0 

DO  70  M=IA(NODE_B) , IA(NODE_A) 
IY=( SLOPE*M+YINT) 
DO  50  L=JA(NODE_B), IY 
IT0T1=IT0T1+IMAGE(M,L) 
IC0UNT1=IC0UNT1+1 
50        CONTINUE 

IF( M.  LT.  IA( NODE_A) )THEN 
DO  60  L=IY+1, JA(NODE_C) 
ITOT2=ITOT2+IMAGE( M, L) 
I COUNT2  =  I COUNT2  + 1 
60         CONTINUE 

END  IF 
70       CONTINUE 

Ll=(N-( IR-1) )*2 
IGREY1=IT0T1/IC0UNT1 
IGREY2=ITOT2/ICOUNT2 
NODE_D=NODE_C+l 

WRITE(20,REC=L1-1)NODE_A,NODE_B,NODE_C, IGREY1 
WRITE(20,REC=L1)NODE_B,NODE_D,NODE_C, IGREY2 
80      CONTINUE 
90     CONTINUE 
RETURN 
END 
C 

Q  ********************************************************** 

SUBROUTINE  EXT_ORIEN( Mil , M12 , M13 , M2 1 , M22 , M23 , 

M31,M32 ,M33 , IENDM, IENDN) 
C 

C        THIS  ROUTINE  DETERMINES  THE  M  MATRIX  PARAMETERS 
C        FOR  ROTATION  OF  THE  IMAGE  PLANE  TO  DESIRED 
C        LOCATION  FOR  VIEWING  IN  THE  SYNTHESIZED  IMAGE. 
C 

IMPLICIT  DOUBLE  PRECISION  (A-Z) 
REAL  MAGN_X/MAGN_Y,MAGN_Z/X/Y,Z 
REAL  X_CORD( 2500 ) , Y_CORD( 2500 ) , Z_CORD( 2500 ) 
REAL  X_VECX, X.VECY, X_VECZ, Y_VECX, Y_VECY, Y_VECZ 
REAL  Z_VECX , Z_VECY , Z_VECZ 
INTEGER  IENDM, IENDN, I TOT, IP, IR 
C 

ITOT=IENDM*IENDN 
DO  10  IR=l,ITOT 
READ( 3,5,END=20)X,Y,Z 
X_CORD(  IR)=X 
Y_C0RD(  IR)=Y 
Z_CORD(  IR)=Z 
5       FORMAT(3(  1X,F17.  7) ) 
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10     CONTINUE 
20     REWIND(3) 

Y_VECX=-( X_CORD(  1 ) ) 

Y_VECY=-( Y_CORD( 1 ) ) 

Y_VECZ=-( Z_CORD(  1 ) ) 

IP=IT0T-IENDM+1 

Z_VECX=X_C0RD( 1)-X_C0RD( IP) 

Z_VECY=Y_C0RD( 1)-Y_C0RD( IP) 

Z_VECZ=Z_C0RD( 1)-Z_C0RD( IP) 
C        USE  THE  CROSS  PRODUCT  OF  Y  CROSS  Z  TO  OBTAIN 
C        THE  X  VECTOR. 

X_VECX=( ( Y_VECY*Z_VECZ)-( Y_VECZ*Z_VECY) ) 

X_VECY=( ( Y_VECZ*Z_VECX)-( Y_VECX*Z_VECZ ) ) 

X_VECZ=(  (  Y_VECX*Z_VECY)-( Y_VECY*Z_VECX) ) 

MAGN_Z=SQRT( ( Z_VECX**2 ) +( Z_VECY**2 ) +( Z_VECZ**2) ) 

MAGN_X=SQRT( ( X_VECX**2 )+( X_VECY**2 )+( X_VECZ**2 ) ) 

MAGN_Y=SQRT(  (  Y_VECX**2 )  +  ( Y_VECY**2 )+( Y_VECZ**2 ) ) 

Ml 1=X_VECX/MAGN_X 

M12=X_VECY/MAGN_X 

M13=X_VECZ/MAGN_X 

M2 1=Y_VECX/MAGN_Y 

M2  2  =Y_VEC Y/MAGN_Y 

M2  3  =Y_VECZ/MAGN_Y 

M3 1=Z_VECX/MAGN_Z 

M3  2 =Z_VECY/MAGN_Z 

M3  3  =Z_VECZ/MAGN_Z 

RETURN 

END 
C 

SUBROUTINE  AFFIN( Al , A2 , Bl , B2 , CI , C2 ) 
C 

C        THIS  SUBROUTINE  ASSIGNS  OR  CALCULATES  THE 
C        COEFFICIENTS  TO  BE  UTILIZED  IN  THE  TRANSFORM  FROM 
C        IMAGE  COORDINATES  TO  SCREEN  COORDINATES. 
C 

IMPLICIT  DOUBLE  PRECISION( A-Z ) 

REAL  Al ,A2, 61,62, CI , C2 ,XIMA_MAX, YIMA_MAX, I_MAX, J_MAX 

DATA  I_MAX,J_MAX/512. 0,512. 0/ 
C 

XIMA_MAX=1600.  0 

YIMA_MAX=1600.  0 

C1=XIMA_MAX 

C2=0. 0 

A2=0. 0 

B1=0. 0 

Al=-XIMA_MAX/( I_MAX*1. 0) 

B2=YIMA_MAX/( J_MAX*1. 0) 

RETURN 

END 
C 

SUBROUTINE  OBS_LOC( XI , Yl , Zl , FOCUS ) 
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c 

C        THIS  SUBROUTINE  CALCULATES  THE  NEW  OBSERVER  X1,Y1,Z1 
C        LOCATION  FROM  DESIRED  LAT.  AND  LONG.  INPUTS  AS  WELL 
C        AS  PROVIDE  THE  FOCAL  LENGTH. 
C 

IMPLICIT  DOUBLE  PRECISION  (A-Z) 

REAL  LATD , LATM , LATS , LOND , LONM , LONS , HE I GHT 

REAL  XI, Y1,Z1,X,Y,Z, FOCUS 
C 

LATD=37. 0 

WRITE(  6,*) *  INPUT  OBSERVER  LATITUDE  IN-MINUTES( REAL) :  ' 

READ(  5,5) LATM 

WRITE(6,*)'  -SECONDS (REAL): ' 

READ( 5, 5) LATS 

LOND=-122.  0 

WRITE( 6,*) ' INPUT  OBSERVER  LONGITUDE  IN-MINUTES( REAL) :  ' 

READ( 5, 5) LONM 

WRITE(6,*)'  -SECONDS (REAL):  ' 

READ( 5, 5) LONS 
5      FORMAT( F5. 1) 

WRITE( 6,*) ' INPUT  OBSERVER  HEIGHT-METERS( REAL) :  ' 

READ(  5, 10) HEIGHT 
10     FORMAT( F6. 1) 

FOCUS=0. 015 

CALL  DMS2XYZ( LATD, LATM, LATS, LOND, LONM, LONS, HEIGHT, 
X,Y,Z) 

X1=X 

Y1=Y 

Z1=Z 

RETURN 

END 
C 

Q      ********************************************************** 

SUBROUTINE  NEW_I J(  XI , Yl , Zl , FOCUS, Al , A2 , Bl , B2 , CI , C2 , 

IENDM, IENDN, IA, JA) 
C 

C        THIS  SUBROUTINE  COMPUTES  THE  NEW  IA  AND  JA  SCREEN 
C        COORDINATES  FROM  THE  GIVEN  OBSERVER  LOCATION. 
C 

IMPLICIT  DOUBLE  PRECISION  (A-Z) 

REAL  XI , Yl , Zl , FOCUS , Mil , M12 , M13 , M21 , M22 , M23 , M3 1 
REAL  XO , YO , X, Y, Z , XIMA , YIMA, Al , A2 , Bl , B2 , CI , C2 
REAL  M32,M33,DENOM 

INTEGER  IENDM, IENDN, ITOT, IR, IA( 2500 ) , JA(  2500) , I , J 
DATA  X0,Y0/0.  0,0.  0/ 
C 

I TOT= IENDN* IENDM 
DO  10  IR=1,2500 
IA( IR)=0 
JA( IR)=0 
10     CONTINUE 

CALL  EXT_ORIEN( Mil , M12 , M13 , M2 1 , M22 , M23 , M3 1 , M32 , M33 , 

IENDM, IENDN) 
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DO  20  IR=1, ITOT 
READ( 3,15, END=30 )X, Y, Z 
15      FORMAT(3(  1X,F17.  7) ) 

DEN0M=M31*(X-X1)+M32*(Y-Y1)+M33*(Z-Z1) 
XIMA=XO-FOCUS*(M11*(X-X1)+M12*( Y-Yl ) +M13*(  Z-Zl ) )/ 

DENCM 
YIMA=YO-FOCUS*(M21*(X-X1)+M22*( Y-Y1)+M23*( Z-Zl) )/ 

DENOM 
XIMA=XIMA* 1000000.  0 
YIMA=YIMA*1000000.  0 

CALL  XY2IJ(XIMA,YIMA,  I , J, Al, A2 , Bl , B2 , CI , C2 ) 
IA( IR)=I 
JA( IR)=J 
20     CONTINUE 
30     REWIND(3) 
RETURN 
END 
C 

Q  *************************************************** 

SUBROUTINE  N0DE_DPTH( XI , Yl, Zl , DEPTH, IENDM, IENDN) 
C 

C        THIS  SUBROUTINE  CALCULATES  THE  DISTANCE  OF  THE 
C        TRANSLATED  ELEVATION  POINTS  TO  THE  NEW  OBSERVER 
C        LOCATION. 
C 

IMPLICIT  DOUBLE  PRECISION( A-Z ) 

REAL  X1,Y1,Z1,DEPTH(2500) 

INTEGER  IENDM, IENDN, IR, ITOT 
C 

I T0T= IENDM* IENDN 

DO  10  IR=1,IT0T 
READ( 3,5,END=20)X,Y,Z 
5       F0RMAT( 3( 1X,F17.  7) ) 

DEPTH( IR)=SQRT( ( (Xl-X)**2)+( ( Yl-Y)**2)+(  (Zl-Z)**2) ) 
10     CONTINUE 
20     REWIND(3) 

RETURN 

END 
C 

Q   ********************************************************** 

SUBROUTINE  FILL( IA, JA, DEPTH, IMAGE, IENDM, IENDN) 
C 

C        THIS  SUBROUTINE  DETERMINES  THE  HIDDEN  SURFACES  AND 
C        CONSTRUCTS  THE  TRANSLATED  IMAGE.  A  Z_BUFFER  IS  USED 
C       TO  HOLD  THE  COMPUTED  DEPTHS  FROM  THE  OBSERVER  TO  THE 
C       GEOGRAPHIC  POSITION.  A  PLANE  EQUATION  IS  CONSTRUCTED 
C        FROM  THREE  ELEV.  POINTS.  THIS  EQUATION  IS  USED  TO 
C        DETERMINE  THE  DEPTH  OF  ALL  POINTS  WITHIN  THE  PLANE. 
C 

IMPLICIT  DOUBLE  PRECISION( A-Z) 

BYTE  IMAGE(512,512) 

REAL  DEPTH( 2500 ) , Z_BUFF( 512 , 512 ) , XI , X2 , X3 , Yl , Y2 , Y3 

REAL  Zl , Z2 , Z3 , Z_DPTH, Cll , C12 , C13 , C21 , C22 , C23 , C3 1 
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REAL  032,033, DET , A_C0EF , B_C0EF , C_C0EF 

INTEGER  IGREY/N0DE_A,N0DE_B,N0DE_C, I_MIN, I_MAX, J_MIN 
INTEGER  J_MAX,FRAME(512,512), IA( 2500 ), JA( 2500 ) , IENDM 
INTEGER  IENDN, IR, I, J,K,L, IPLANES 
C 

C        FIRST  DETERMINE  THE  NUMBER  OF  PLANES  AND  INITIALIZE 
C        THE  Z_BUFFER  AND  IMAGE  TO  0. 
C 

IPLANES=( ( IENDM-1)*2)*( IENDN-1) 
DO  20  K=l,512 
DO  10  L=l,512 
IMAGE(L,K)=0 
Z_BUFF( L,K)=0 
10      CONTINUE 
20     CONTINUE 
C 

C        DETERMINE  THE  COEFFICIENTS  OF  THE  PLANE  EQUATION 
C        FROM  THE  THREE  ELEVATION  POINTS. 
C 

DO  70  IR=1, IPLANES 
READ ( 2  0 , REC= I R ) NODE_A , N0DE_B , NODE_C , I GREY 
Xl=IA(NODE_A) 
Yl=JA(NODE_A) 
Zl=DEPTH(NODE_A) 
X2=IA(NODE_B) 
Y2=JA(NODE_B) 
Z2=DEPTH(NODE_B) 
X3=IA(NODE_C) 
Y3=JA(NODE_C) 
*   Z3=DEPTH(NODE_C) 
C        DETERMINE  THE  COFACTOR  ELEMENTS 
Cll=( ( Y2*Z3)-(Y3*Z2) ) 
C12  =  -(  (X2*Z3)-(X3*Z2)) 
C13=(  (X2*Y3)-(X3*Y2) ) 
C21=-( ( Yl*Z3)-( Y3*Z1)) 
C22=( (X1*Z3)-(X3*Z1)) 
C23=-(  (X1*Y3)-(X3*Y1) ) 
C31=( ( Y1*Z2)-(Y2*Z1) ) 
C32=-( (X1*Z2)-(X2*Z1) ) 
C33  =  (  (X1*Y2)-(X2*Y1) ) 
C        CALCULATE  THE  DETERMINANT 
C 

DET=( X1*C11 )+( Y1*C12 ) +( Z1*C13 ) 
C 

C       THE  COEFFICIENTS  ARE  DETERMINED  FROM  MULTIPLYING 
C       THE  reciprocal  of  the  determenant  with  the 
C       transpose  of  the  cofactors  called  the  adjoint. 
C 

A_C0EF=-( (C11+C21+C31)/DET) 
B_C0EF=-( ( C12+C22+C32 )/DET) 
C_C0EF=-( (C13+C23+C33)/DET) 
IF(C_C0EF.EQ.  0.0)GOTO  70 
C        DETERMINE  THE  MAXIMUM  AND  MINIMUM  VALUES  TO  TEST 
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C        FOR  THE  FILL  ALGORITHUM. 
C 

I_MAX=MAX( IA( NODE_A) , IA( NODE_B) , IA( NODE_C) ) 
I_MIN=MIN( IA(NODE_A), IA(NODE_B), IA(NODE_C) ) 
J_MAX=MAX( JA( NODE_A) , JA( NODE_B ) , JA( NODE_C ) ) 
J_MIN=MIN(  JA(  NODE_A) , JA( NODE_B ) , JA( NODE_C ) ) 
IF( I_MIN. LT. 1. OR. I_MAX. GT. 512 )GOTO  70 
IF( J_MIN. LT. 1. OR. J_MAX. GT. 512 )GOTO  70 
C 

C        CLEAR  THE  REFERENCE  FRAME  BUFFER  AND  CALL  THE 
C        FRAME  FILL  SUBROUTINE. 
C 

DO  40  L=I_MIN, I_MAX 
DO  30  K=J_MIN, J_MAX 
FRAME( L,K)=0 
30       CONTINUE 
40      CONTINUE 

CALL  FRAME_FIL(NODE_A,NODE_B,NODE_C, FRAME, IA,JA, 

J_MAX,J_MIN) 
DO  60  J=J_MIN, J_MAX 
DO  50  I=I_MIN, I_MAX 
IF(FRAME( I, J). EQ. 1)THEN 

Z_DPTH=-( l+(A_COEF*I)+(B_COEF*J) )/C_COEF 
IF( Z_BUFF( I, J). EQ. 0. 0. OR. Z_DPTH. LT. Z_BUFF( I, J) )THEN 
Z_BUFF( I ,  J ) =Z_DPTH 
IMAGE( I, J)=IGREY 
END  IF 
END  IF 
50       CONTINUE 
60      CONTINUE 
70     CONTINUE 

DO  90  J=l,512 
WRITE(21,REC=J)( IMAGE( I, J), 1=1,512) 
90     CONTINUE 
RETURN 
END 
C 

Q  ********************************************************* 

SUBROUT INE  FRAME_F I L( NODE_A , NODE_B , NODE_C , FRAME , 

IA, JA, J_MAX, J_MIN) 
C 

C        THIS  SUBROUTINE  CONSTRUCTS  AN  EDGE  LIST  FROM  THE 
C        THREE  NODES  PASSED  IN  THE  ROUTINE.  THESE  EDGES 
C        ARE  USED  IN  A  POLYGON  FILL  ROUTINE  USEING  A  FRAME 
C       BUFFER  AND  Y_BUCKET 

IMPLICIT  DOUBLE  PRECISION( A-Z) 

REAL  I_INCPT,DELTA_I,DELTA_J,AEL( 3,4) 

REAL  XNODEl/XNODE2/DX 

INTEGER  NODE_A, NODE_B, NODE_C, FRAME ( 512 , 512 ) , IA( 2500 ) 

INTEGER  JA(2500), J_MAX, J_MIN/NODE(4)/ IS, NODE1 , NODE2 

INTEGER  N_HIGH,N_LOW, IT, IU, J_BUCKET( 512), IR, II, 12 

INTEGER  ICNT,LCNT 
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DO  10  IS=1,512 
J_BUCKET( IS)=0 
10     CONTINUE 

DO  30  IS=1,3 

DO  20  IR=1,4 
AEL(  IS, IR)=0.  0 
20      CONTINUE 
30     CONTINUE 

NODE( l)=NODE_A 
NODE(2)=NODE_B 
NODE( 3)=NODE_C 
NODE( 4)=NODE_A 
DO  40  IS=1,3 
NODEl=NODE(  IS) 
NODE2=NODE(  IS+1) 
IF(  JA(  NODE1 ) .  GE.  JA(  NODE2  )  )THEN 
N_HIGH=N0DE1 
N_LOW=NODE2 
ELSE 
N_HIGH=NODE2 
N_LOW=NODEl 
END  IF 

I_INCPT=IA( N_HIGH) 
DELTA_J=( JA(N_HIGH)-JA(N_LOW) ) 
IF(DELTA_J. EQ.  0. 0)THEN 

GOTO  40 
ELSE 

DELTA_I=-( IA(N_HIGH)-IA(N_LOW) )/DELTA_J 
END  IF 

IT=JA(N_HIGH) 
IU=J_BUCKET(  IT) 
IF( IU. EQ. 0)THEN 

J_BUCKET(  IT)  =  IS 
ELSE 

AEL( IU/4)=IS 
END  IF 

AEL( IS/1)=I_INCPT 
AEL( IS,2)=DELTA_I 
AEL( I S , 3 ) =DELTA_ J 
40     CONTINUE 

IT=J_BUCKET( J_MAX) 
IU=INT(AEL(  IT, 4)) 
XN0DE1=AEL(  IT,1) 
XNODE2=AEL(  IU,1) 
DX=XN0DE1-XN0DE2 
IF(DX.  LE.  0.  0)THEN 
I1=NINT(XN0DE1) 
I2=NINT(XNODE2) 
ELSE 
I1=NINT(XN0DE2) 
I2=NINT(XN0DE1) 
END  IF 
DO  50  IR=I1, 12 
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FRAME( IR, J_MAX)=1 
50     CONTINUE 

AEL( IT/3)=AEL( IT,3)-1. 0 

AEL( IU,3)=AEL( IU,3)-1. 0 
ICNT=J_MAX-1 
LCNT=J_MIN 

DO  70  IS=ICNT/LCNT/-1 
I F( J_BUCKET(  I S ) .  EQ.  0 ) THEN 
XN0DE1=XN0DE1+AEL(  IT, 2) 
XNODE2=XNODE2  +  AEL(  IU,  2  ) 
ELSE  IF(AEL( IT, 3). LE.  0. 0)THEN 
IT=J_BUCKET(  IS) 
XN0DE1=AEL(  IT, 1) 
XNODE2=XNODE2+AEL(  IU,2) 
ELSE 
IU=J_BUCKET(  IS) 
XN0DE1=XN0DE1+AEL(  IT, 2 ) 
XN0DE2=AEL(  IU, 1) 
END  IF 

DX=XNODE 1 -XN0DE2 
IF(DX.  LE.  0.  0)THEN 
I1=NINT(XN0DE1) 
I2=NINT(XN0DE2) 
ELSE 
I1=NINT(XN0DE2) 
I2=NINT(XN0DE1) 
END  IF 

DO  60  IR=I1, 12 
FRAME ( IR, IS)=1 
60      CONTINUE 

AEL(  IT,3)=AEL(  IT,3)-1.  0 
AEL( IU,3)=AEL(  IU,3)-1.  0 
70     CONTINUE 
RETURN 
END 
C 

Q   ********************************************************* 
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